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ABSTRACT 
AnonHlous  transport  caused  by  fluctuations  of  resistive  pressure- 
driven  modes  is  discussed  within  the  framework  of  magnetohydrodynamics 
(MHD).  The  nonlinear  reduced  equations  describing  fluctuations 
localized  near  a  particular  magnetic  field  line  are  derived  for  tokaraak 
and  reversed-field -pinch  (RFP)  plasmas,  taking  into  account  nonzero 
viscosity  and  heat  conductivity.  Based  on  the  reduced  equations,  the 
relationship  between  the  dynamo  effect,  or  the  electric  field  caused  by 
the  fluctuations,  and  the  anomalous  heat  transport  is  obtained.  For  an 
ideally  stable  but  resistively  slightly  unstable  plasma,  the  anomalous 
transport  is  caused  particularly  by  convective  motions.  The  convection 
is  studied  as  bifurcation  from  the  linearly  unstable  equilibrium  and 
the  expression  of  the  anomalous  transport  in  a  tokamak  plasma  is 
obtained  as  a  function  of  the  mean  pressure  gradient  near  the  critical 
point.  In  order  to  evaluate  the  effects  of  the  convection  on  the 
anomalous  transport  under  various  conditions,  the  reduced  equations  are 
also  solved  numerically.  It  is  found  that  the  Nusselt  number,  that  is, 
the  ratio  of  the  total  heat  conductivity  including  the  anomalous  heat 
transport  to  the  classical  collisional  heat  conductivity,  is 
significantly  large  under  some  conditions.  This  partially  explains  the 

large  heat  losses  in  controlled  thermonuclear  fusion  devices. 
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I.   INTRODUCTION 

In  a  magnetically  confined  plasma,  various  types  of  fluctuations 
affect  the  global  behavior  of  the  plasma  through  enhancement  of  the 
transport  of  mass,  heat  and  magnetic  fields.  The  collisions  of  ions 
and  electrons  always  generate  some  diffusion  in  plasmas,  which  is 
called  the  collisional  diffusion.  However  the  transport  caused  by  the 
fluctuations,  which  we  call  the  anomalous  transport,  has  in  many  cases 
similar  or  even  greater  effects  on  the  plasma  than  the  collisional 
diffusion.  The  anomalous  transport,  therefore,  has  been  a  subject  of 
much  interest  in  recent  years. 

In  this  thesis,  we  will  consider  "low-level"  or  "weak" 
fluctuations  or  the  fluctuations  with  small  amplitude  which  vary  on 
faster  time  scales  and  smaller  spacial  scales  than  the  mean  fields  in  a 
plasma.  The  fluctuations  commonly  observed  in  well  confined  plasmas 
are  always  weak,  typically  It  of  magnetic  field  fluctuations  in 
tokamaks  and  ^0%  in  reversed  field  pinches  (RFP's).  The  steady 
convection  or  the  saturated  islands  on  a  rational  magnetic  surface  and 
plasma  turbulence  are  examples  of  the  weak  fluctuations  that  we  are 
particularly  interested  in. 

We  will  exploit  the  weakness  of  the  fluctuations  in  a  plasma  in 
order  to  derive  the  set  of  equations  of  the  fluctuations.  For  general 
turbulence  theories  in  fluid  dynamics,  many  interesting  problems  are 
associated  with  reasonably  strong  fluctuations.  In  such  cases,  it  is 
known  that  the  system  of  equations  describing  fluctuating  quantities 
would  not  be  closed  without  further  approximation  (such  as  the  closure 
method^).  The  model  of  weak  fluctuations  to  be  employed  in  this 
thesis,  however,  enables  us  to  circumvent  these  difficulties.  In  this 
model,  introducing  a  consistent  ordering  of  every  physical  quantity,  we 
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are  able  to  close  the  system  to  the  lowest  order.  For  example,  the 
nonlinear  term  in  Faraday's  law  for  the  perturbed  magnetic  field  B,  is 
7x(v,xB,  -  <v,xb^>)  where  v^  is  the  perturbed  velocity  field.  Since 
the  derivative  of  term  v^^B^  is  much  larger  than  the  derivative  of  its 
average  <v,'<B^>,  where  <  >  denotes  the  average  of  the  fluctuations 
which  smoothes  away  the  large  derivatives  of  the  fluctuating 
quantities,  the  term  Vx<v^xB^>  does  not  appear  in  the  equation  for  B, 
to  the  lowest  order.  This  procedure  significantly  simplifies  the 
problem. 

One  of  the  most  intriguing  problems  which  is  believed  to  be 
related  to  the  weak  fluctuations  in  a  magnetically  confined  plasma  is 
the  spontaneous  reversal  of  the  toroidal  magnetic  field  in  RFP 
experiments.^  Taylor  suggested  that  a  slightly  resistive  plasma 
minimizes  its  energy  through  a  magnetic  reconnection  process,  subject 
to  the  constraint  that  the  total  magnetic  helicity  of  the  plasma  be 
conserved.  Although  he  did  not  specify  any  particular  dynamical 
process  of  the  plasma  relaxation,  the  predicted  minimum  energy  state 
agrees  well  with  the  experimentally  observed  equilibrium  state  of  an 
RFP  plasma.  Recently  it  has  been  demonstrated^ '^ '^  that  the  dynamical 
description  of  a  plasma  based  on  the  nonlinear  resistive  MHD  equations 
with  the  presence  of  weak  fluid  fluctuations  successfully  describes  the 
relaxation  of  a  plasma  to  a  state  similar  to  the  one  predicted  by 
Taylor. 

Another  interesting  phenomenon  associated  with  the  weak 
fluctuations  is  the  anomalous  heat  transport.  For  instance,  turbulence 
in  an  RFP  and  steady  convection  on  some  rational  magnetic  surfaces  in  a 
tokamak  enhance  the  heat  conduction  across  the  magnetic  flux  surfaces 
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and  deteriorate  energy  confinement  of  the  plasmas.  Thus  it  is 
important  to  understand  how  the  anomalous  heat  transport  behaves  for  a 
given  set  of  macroscopic  profiles  of  the  plasma,  such  as  the  pressure 
profile,  in  order  to  determine  the  energy  confinement  properties  of  the 
plasma.  Determination  of  the  dependence  of  the  anomalous  heat 
transport  on  such  macroscopic  conditions  is  one  of  the  main  goals  of 
this  thesis. 

The  class  of  the  weak  fluctuations  with  which  we  are  particularly 
concerned  here  is  the  one  generated  by  the  resistive  g-mode  or,  more 

Q 

precisely,  the  resistive  fast  interchange  mode,  which  is  the 
instability  caused  by  pressure  or  temperature  gradients  acting  against 
the  curvature  c"  -.e  magnetic  field  lines  in  a  plasma  with  finite 
resistivity,  .  .gous  to  the  thermal  instability  causing  the  Benard 
convection  in  a  fluid  heated  from  below  in  a  gravity  field.  These 
resistive  pressure  driven  modes  are  expected  to  be  present  in  tokamaks 
and  RFPs.  Recently,  Hameiri'  discussed  the  turbulent  heat  conduction 
without  specifying  any  particular  modes  and  derived  some  of  its  general 
properties.  In  this  thesis  we  will  consider  the  anomalous  heat 
transport  specifically  caused  by  the  resistive  g-mode  fluctuations  and 
we  will  determine  a  more  precise  characterization  of  the  transport. 

For  the  resistive  g-mode,  the  free  energy  source  is  the  mean 
pressure  (or  temperature)  gradient,  and  the  energy  sink  is  the 
collisional  diffusion.  Therefore  when  the  mean  pressure  gradient 
exceeds  a  critical  value,  modes  localized  on  a  rational  surface  begin 
to  grow.  These  modes,  however,  eventually  saturate  with  finite 
amplitude  if  the  mean  pressure  gradient  is  not  too  large  compared  to 
the  critical  value,  and  we  have  steady  convection  within  a  boundary 
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layer  which  is  analogous  to  the  Benard  convection  in  fluid  dynamics. 
This  resemblance  of  the  resistive  g-mode  to  the  Benard  convection  was 
first  pointed  out  and  studied  by  Dagazian^  based  on  a  quasi-linear 
theory  °  with  the  shape  assumption  or  the  assumption  that  in  a  slightly 
nonlinear  regime  the  mode  retains  its  linear  shape  with  a  finite 
amplitude.  We  will  deal  here  with  the  similar  problem  of  the 
convection  arising  from  the  resistive  g-mode  using,  however,  a 
different  method.  Our  approach  is  to  use  nonlinear  bifurcation 
analysis  to  derive  the  dependence  of  the  anomalous  heat  transport  on 
the  mean  pressure  gradient  analytically.  In  this  method,  the  set  of 
the  nonlinear  equations  of  the  fluctuations  is  expanded  in  terms  of  the 
small  amplitude,  and  the  complete  algorithm  to  determine  all  the  higher 
order  terms  is  obtained.  Therefore  we  do  not  need  any  further 
assumptions  such  as  the  shape  assumption.  We  will  also  show  the 
results  of  numerical  computations  of  the  nonlinear  equations  of  the 
fluctuations,  which  determine  the  dependence  of  the  anomalous  heat 
transport  on  the  mean  quantities  under  various  conditions. 

This  thesis  consists  of  nine  sections.  Section  II  reviews  and 
extends  the  set  of  equations  governing  the  mean  quantities  obtained  by 
Hameiri.'  In  Section  III,  we  derive  the  nonlinear  reduced  equations 
describing  the  fluctuations  due  to  the  resistive  g-modes  with  finite 
viscosity  and  heat  conductivity,  for  both  an  RFP  and  a  tokamak. 
Sections  IV  and  V  precsent  general  discussions  on  the  anomalous 
transport  in  an  RFP  and  a  tokamak,  respectively.  Since  the  set  of 
reduced  equations  for  a  tokamak  is  simpler  and  easier  to  treat,  Section 
VI  discusses  its  nonlinear  stability,  from  which  we  see  the  similarity 
of  the  resistive  g-mode  and  the  Benard  convection.   Section  VII 
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contains  the  bifurcation  analysis  of  steady  convection  of  a  tokamak 
plasma.  Here  the  dependence  of  the  anomalous  transport  due  to  the 
convection  on  the  mean  pressure  gradient  is  derived  near  the  critical 
pressure  gradient.  Section  VIII  summarizes  the  aforementioned 
properties  of  the  anomalous  transport  with  some  quantitative  results 
obtained  from  the  numerical  simulations.  Finally,  Section  IX  contains 
discussions  and  conclusions. 
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II.      EQUATIONS   OF   THE  MEAN  MOTION 

We     begin     by     reviewing      the     derivation     of  the  set  of  equations 
describing   the  evolution  of   the  mean  quantities.        Here     we     will     also 

generalize     the     earlier     work^     by     taking      into     account   the  parallel 

1  2 
viscosity       in  addition  to  the  other  diffusion  coefficients. 

We  start  from   the   following  resistive  MHD  equations: 


p(_.  +   y«V]v   -   -Vp  +   JxB   -   7-n  (2-1a) 


3B 

-=.  *  Vx(nJ   -  vxB)   -  0  (2-1b) 

o  t  —  -     — 

1^  *   V.(pv)    -   0  (2-lc) 

_£  +    v'Vp   ■•■   YpV'V   -    (Y-1)(nJ'J    +   V«(<7T)    -   n:Vv)    ,  (2-1d) 

ot~  ~  ~~  »~ 


where 


P  =   pT  (2-1  e) 

n   -  -3y.  .(bb  -  1  DA   -   y  o  (2-lf) 

A    -   b-[(b-V)v]    -  1  V-y  (2-1g) 

b   -   B/|b|  (2-1h) 


In  Eqs.  (2-1),  symbols  have  their  usual  meanings:  p,  p,  T  and  v 
are  the  mass  density,  pressure,  temperature  and  velocity,  respectively. 
B  is  the  magnetic  field,  satisfying  V«B  -  0,  J  -  7x8  is  the  current 
density,  Y  is  the  ratio  of  the  specific  heats,  n  is  the  resistivity 
tensor,    and  <  is   the  heat     conductivity     tensor.        Both     n     and     <     are 
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assumed  to  be  positive  definite.  Equation  (2-1 f)  is  the  definition  of 
the  stress  tensor  n,  where  I  denotes  the  unit  tensor  and  o  is  the  rate- 
of-strain  tensor  defined  by 


3v^   3v.   2 
iJ   3xj   3x^   3  ^J  - 


for  the  (i,j)  component  in  Cartesian  coordinates.   Here  6j^j   indicates 

Kronecker's  delta,   that  is,   6^-1   if  i  -  j  and  6^j  -  0  if  i  *  j. 

V-n  =  I   311.  ./3x.  is  the  divergence  of  the  tensor  n  and  IliVy  is  the 

i 

contraction  defined  as  n:Vv  =  11   I1^.3v^/3x..   It  follows  that 

i  j 


n:Vv  =  -3v      \^  -  ^  tr    0^ 


where  tr  denotes  the  trace  of  the  tensor. 

We  assume  that  the  plasma  is  confined  in  either  a  toroidal  or  a 
cylindrical  vessel  with  a  perfectly  conducting  wall.  Therefore  we 
impose  the  boundary  conditions 

n«B  -  0 

V  -  0 
and 

Exn  -  0  , 

where  n  denotes  the  unit  vector  normal  to  the  boundary,  and  E  denotes 
the  electric  field,  defined  by 
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E  -  nJ  -  vxB  . 

In  case  of  a  cylindrical  vessel,   we  identify  the  two  ends  of  the 
cylinder  and  require  periodicity  of  all  physical  variables. 

Equations  (2-1  a)  and  (2-1  d)  may  be  written  in  the  equivalent  form: 


—  (pv)  +  V«(pvv)  -  -Vp  +  VxB  -  7-n  (2-2) 

9t   -        --  -  - 


^   (ps)  +  V'(psv)  =  -   nJ'J  +  V«[l  <Vt1  +  —  k7T'7T 


3t  ^  -    J  -  _     'T        j2 


*|u.  X2  *±M 
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I 


tr  0^  .  (2-3) 


Here  s  denotes  the  specific  entropy  defined  as 


s  -  J-  log  (p/p^)  .  (2-iJ) 


For  convenience,  we  work  with  variables  measured  in  two  natural 
units,  in  which  distances  are  measured  in  units  of  the  plasma  minor 
radius  a,  times  are  measured  in  units  of  the  Alfven  transit  time 
a  /  p/B,  where  B  and  p  are  typical  values  of  the  magnitude  of  the 
magnetic  field  and  the  mass  density  of  the  plasma,  and  the  magnetic 
field  and  mass  densities  are  measured  relative  to  B  and  p, 
respectively.  Other  variables  are  then  measured  accordingly. 

We  now  consider  the  fluctuations  in  the  plasma  described  by  Eqs. 
(2-1).  Since  we  are  interested  in  evolution  of  "mean"  physical 
quantities  under  some  influence  of  turbulence,  we  have  to  clarify  what 
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we  mean  by  a   "mean"  or   "averaged"  value  of  each  physical  quantity.     Let 
U3     denote     <f>      the     averaged  value  of  a  physical  quantity   f.      Then  we 

require  that   the  averaging  operator  <   >   satisfy   the  so-called     Reynolds 

2 
conditions  : 

<f+g>  -  <f>  +  <g>  ,  (2-5a) 

<af>  -  a<f>  ,  (2-5b) 

<a>  =  a  ,  (2-5c) 

<">  .  ^  .    (,  <  1  <  3)  ,    <|f>  ■  ^  .         (2-5<l) 

9x^      dx^  ~        ~  dt      dt 

<<f>g>  =  <f><g>  ,  (2-5e) 

where  f  and  g  are  some  physical  quantities  and  a  is  a  constant.  Since 
we  consider  low-level  fluctuations  varying  rapidly  in  space  and  time 
superimposed  on  a  mean  motion,  which  varies  slowly  in  space  and  time, 
we  may  think  of  <  >  as  either  an  ensemble  average  or  an  average  over 
the  small  space-  and  fast  time-scales  of  the  fluctuations.  However, 
the  following  discussion  in  this  section  does  not  depend  on  the 
particular  choice  of  the  averaging  operation  <  >  so  long  as  it 
satisfies  conditions  (2-5).  Therefore  we  do  not  specify  the  actual 
form  of  the  averaging  here  and  will  wait  until  further  assumption  is 
made  on  the  fluctuations  so  that  we  can  easily  specify  the  most 
appropriate  form  of  averaging. 

Once  an  averaging  operator  <  >  satisfying  conditions  (2-5)  is 
defined,  we  can  write  every  physical  quantity  as  the  sum  of  the  mean 
part,  denoted  by  subscript  o  and  the  fluctuating  part,  denoted  by 
subscript  1 .   This  decomposition  is  unique.  For  example,  the  magnetic 
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field  §  -  §0  ■*■  §!•  where  Bq  -  <B>  and  B^  -  B  -  <B>.  It  is  easy  to  show 
from  properties  (2-5)  that  <B^>  -  0. 

For  simplicity,  we  take  n  to  be  a  small  scalar  quantity.  We  now 
derive  the  dynamical  equations  of  the  mean  quantities.  Since  we  are 
dealing  with  low-level  fluid  fluctuations,  we  use  the  following 
assumptions:  (1)  mean  quantities,  taken  to  be  0(1)  except  for  the  mean 
f"low  Vq,  evolve  on  a  resistive  time  scale  and  vary  on  the  minor  radius 
length  scale,  i.e.,  alv^l  =  0(n),  af^  ■—  '   0(n)  and  a|V|  «  0(1);  (2)  the 

"  o  t 

fluctuations  may  vary  more  rapidly  in  space  and  time  than  the  mean 
quantities,   i.e.,   a'^  ■—  >>  0(n),  a|V|  >>  0(1);  (3)  the  energy  of  the 

o  t 

fluctuations  is  considerably  smaller  than  the  mean  magnetic  energy, 
i.e.,  |bJ2,  PoI^iI^.  Pi-  Pq"^!'  '^0^1'  '^^  l§ol^-  We  note,  therefore, 
that  the  fluctuating  current  J^  •  VxB^  is  of  the  order  of  J^,  or 
^l^ll^l§ol  "  0(1);  (^)  the  perpendicular  viscosity  u  is  taken  to  be 
0(ri)  while  the  parallel  viscosity  y  is  taken  to  be  0(1).  Here  we 
also  only  consider  the  perpendicular  component  <  and  the  parallel 
component  <  of  the  conductivity  tensor  <:  it  is  assumed  that  <  and 
<  are  0(ri)  and  0(1),  respectively.  In  general,  these  diffusion 
coefficients  n,  u  and  <  depend  on  such  physical  variables  as  p ,  T  and  B 
and,  therefore,  they  also  fluctuate.  However,  from  these  ordering 
assumptions,  it  is  easy  to  check  that  the  fluctuating  parts  of  the 
diffusion  coefficients  are  so  small  that  they  enter  neither  the 
equation  of  the  mean  motion  nor  the  equation  of  the  fluctuations. 

Under  these  ordering  assumptions,  we  obtain  the  set  of  equations 
for  the  mean  quantities  by  taking  the  average  of  the  MHD  equations. 
Averaging  {2-2).  (2-1b),  (2-1c),  (2-3),  (2-1e)  and  {2-H)  yields, 
respectively. 


-11- 


3B 

-I^  +  Vx(nJ„  -  v^xB^  -  <v.xBi>)  -  0  (2-6b) 

3t       ~°   ~°  ~°    "'  ~ 

3p 

— 2  +  V'(pyQ  +  <PiYi>)  ■  0  (2-6c) 

A  (p^Sq)*  7.(Po3qVq  *  Po<3iVi>  *  3^<PtVt>) 

=  J-  (nJ^-Jo  *  <nJi-Ji>)  *  V.[J-  <VtJ  *  J-  (<VT^.VTo  *  «VT,.7Ti>) 
^o  o         o 

+  J_  u  <A?>  +  —  u  tr  <o?>  (2-6d) 


To   II   ^     2T^  1 


p  -  p  T  (2-6e) 


^o  =  1°8  7^  (Pq/pJ^  •  ^2-6f) 


Here  Eqs.    (2-6b)  and  (2-6c)  are  exact  while  Eqs.   (2-6a),  (2-6e)  and 
(2-6f)  holds  to  0(1),  Eq,   (2-6d)  to  0(n). 

There  is  an  important  relationship  among  the  averaged  quantities: 


<Yi»Bi>-Jo  *  TTTTyr  'Pi-i'*'^o  -  ^o  <P^l^"'^^o 


-   -    (<^Jr^1>   *  —   <<VT,.VT,>   *    3li      <Af>   *  -f  tr   <of>]    *  f        (2-7) 
"     ~  ^0  I  I 


where 
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f   -   -V-<(p,    +   BqB,)v,    -    (v^.B,)Bq   ♦  I  Pov2.v>    .  (2-8) 


This  may  be  derived  as  follows:    the  energy  equation  of     the     system     is 
found  by  combining  Eqs.      (2-1) 

9     r 1    n2         1         2  P    1 

—   I_  B      +  -   pv      +  — !— J 

9t      2    -  2      -  Y-1 

+   V.[(nJ   -   vxB)xB   +  -  pv^v   +  pv   -   kVT   -   IIv]    -    0    . 

-        -    -      -        2      -    -        Y-1       - 

Averaging  this  equation,  we         obtain       the       expression       for 

9(Bq/2  +   pQ/(Y-1)]/3t   to  order  of  n.      We  also  form     another     expression 
for       3[Bq/2  +   pQ/(Y-1)j/3t        from       the       averaged       equations     (2-6). 
Subtracting    these  two  expressions  yields  Eq.      (2-7). 
By  using   the   following  definitions 


c   -   <v,xB.>   -  -L  <p.v^>xB-  (2-9) 

-1-1  p  1-1        -0 


Ho  •  ^0  *7^  <PlYi>    .  (2-10) 

^o 


Equation  (2-7)  may  be  written  as 


£•^0  *  Po^^l^^l^'^'^o 


-  -  {<n'l^-J^>   ♦  -L  <k7Ti'VT,>  +  3y   <Af>  +  -^  tr  <af>}  *  f  ,   (2-11) 
~  ~     ^0  II 
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where  the  relation  s-^    -  P^ /Pq  "  ^Pi^Pq  "^^  used.   Writing  Eq.    (2-6d) 
in  terms  of  p^  and  using  Eq,   (2-11),  we  rewrite  Eqs.   (2-6): 

3B 

-Zl  *   7«(nJo  -  £  -  UqxBo)  -  0  (2-12a) 

3Po 

=  (ir-l)[nJo-Jo  -  £-Jo  +  V.(<VTq)  -  7-(Pq  <3iV,>)  *  fj       (2-12C) 

From  these  equations,  it  is  clear  that  e  and  -pQ<3^y^>  are  the  extra 
electric  field  and  the  extra  heat  flux  caused  by  the  fluctuations, 
respectively.  We  therefore  call  e  the  anomalous  electric  field  and 
-p  <s,v^>  the  anomalous  heat  flux.  We  also  call  the  effect  of  £  the 
dynamo  effect  since  it  is  related  to  the  generation  of  the  magnetic 
field. '"'^ 

Together  with  Eq.  (2-6a),  Eqs.  (2-12)  -re  the  closed  set  of 
equations  to  evolve  the  mean  quantities  if  the  anomalous  electric  field 
e,  the  anomalous  heat  flux  -Po*"SiV^>  and  the  function  f  are  given.  We 
note  that  f  vanishes  for  localized  fluctuations  as  will  be  discussed. 
We  also  note  that  the  quantity  Uq  is  determined,  as  in  the  Grad-Hogan 
model  for  diffusion,  by  the  requirement  that  all  mean  quantities 
evolve  through  a  sequence  of  magnetostatic  equilibrium  states 
satisfying  Eq.  (2-6a).  The  dependence  of  the  anomalous  electric  field 
e  on  the  mean  profiles  was  determined  by  Bhattacharjee  and  Hameiri,  ' 
for     the     case  of  tearing-mode-induced  turbulence.      One  of  the  goals  of 
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this   thesis  is  to  determine   the  anomalous  heat  flux     from     the     set     of 
equations     of     the  resistive  g-mode  fluctuations,   which  will  be  derived 
in  the  next  section. 

We  now  briefly  comment  on  some  properties  of  the  anomalous 
electric  field  e  and  the  anomalous  heat  flux  -Pq<s^v^>.  By  integrating 
Eq.      (2-11)  over  the   total  volume  of   the  plasma,    we  obtain 


^^'io  *   Pq^^I^I^'^'^'o^'^^^ 


=   -J    [<nJT'Ji>   +  —   <<VT,'VT,>   +    3u      <xf>   +  -1  tr   <a^>]6^    x    ,      (2-13) 
~     ~  0  II 

since  the  volume  integral  of  f  vanishes.  Equation  (2-13)  gives  the 
relationship  of  the  integrals  of  e  and  <3iVi>  over  the  total  plasma 
volume.  Moreover,  in  certain  limiting  cases,  the  term  e'^q  appearing 
in  Eqs.  (2-11)  and  (2-1 2e),  is  shown  to  vanish.  In  fact,  Hameiri 
proved'  that  e'Jq  -  0  if  the  resistivity  n  -  0.  We  will  also  show  in 
Section  V  that  the  term  t'J^  vanishes  if  the  pressure  gradient  VPq  is 
taken  to  be  vanishingly  small.  In  either  of  these  cases,  we  have,  from 
Eq.      (2-13), 


J    pQ<s^y,>.VTQd3x   <  0  (2-li4) 


This  inequality  implies  that  the  anomalous  heat  flux  tends  to  transport 
heat  in  the  direction  opposite  to  the  mean  temperature  gradient  VT^. 
This  is  similar  to  the  collisional  heat  conductivity.  As  a  special 
case,  if  the  fluctuations  are  highly  localized  and  if  we  may  take  the 
integral     in     a     small     volume,      e.g.,      in     the     vicinity  of  a  rational 
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surface,      we     expect     that     Pq'^^i^i'* '^"^o   *"  °       holds       locally.  This 

consideration  suggests  that,  in  this  case,  the  anomalous  heat  flux 
-Pq<3^v^>  has  a  form  of  K^VT^,  where  K^  is  a  non-negative  function 
which  may  depend  on  mean  quantities,    including  VTq. 

In  the  following  sections,  we  will  discuss  the  fluctuations  due  to 
the  resistive  g-mode  and  characterize  the  anomalous  heat  flux  in  more 
detail. 
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III.   REDUCED  EQUATIONS  OF  THE  FLUCTUATIONS 

We  will  now  derive  the  equations  describing  the  dynamics  of  the 
fluctuating  quantities  arising  from  the  resistive  g-mode.  This  is  an 
extension  of  the  reduced  equations  obtained  in  Ref.   1. 

We  now  assume  that  the  mode  is  localized  along  a  mean  magnetic 
field  line  so  that  the  perpendicular  derivative  of  a  fluctuating 
quantity,  i.e.,  the  derivative  across  the  field  line  is  of  order  1/6, 
where  6  is  a  small  parameter  measuring  the  localization  of  the  mode, 
while  the  parallel  derivative  of  a  fluctuating  quantity,  i.e.,  the 
derivative  along  the  field  line,  is  of  order  1.  Assuming  that  the 
resistivity  n  is  a  scalar  constant  for  simplicity,  we  take  6  to  be 
0(/n).  The  time  derivative  of  a  fluctuating  quantity  is  assumed  to  be 
0(1),  while  the  time  derivative  of  a  mean  quantity  is  0(n)  as  mentioned 
before.  All  the  mean  quantities  are  at  most  of  order  1  except  for  v^ 
which  is  0(6  ),  and  the  fluctuating  quantities  v^ ,  B^ ,  p^  and  T^  are 
assumed  to  be  0(6).  These  orderings  are  consistent  with  the  assumption 
we  made  in  Section  II  and  agree  with  the  scaling  of  the  resistive  fast 

Q 

interchange  mode  in  the  linear  stability  theory. 

We  now  describe  the  nonlinear  mode  equations  according  to  these 
orderings.  For  simplicity  we  take  the  parallel  viscosity  coefficient 
M  to  be  smaller  that  0(1)  so  that  it  does  not  eneter  the  equation  of 
the  fluctuations.  The  more  complete  version  of  the  mode  equations 
including  the  parallel  viscosity  coefficient  y  is  derived  in  Appendix 
A.  The  perpendicular  viscosity  is  taken  to  be  0(6  ),  as  before. 
Subtracting  Eq.   (2-6a)  from  Eq.   (2-1  a),  we  obtain  to  0(6) 


po^^  *  Yr^hi  ■  -^Pi  *  io"?!  *  ^r5o  *  ^r§i  *  Yi-^  *       ^^"^^ 
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Similarly,    subtracting  Eq.      (2-6b)    from  Eq.      (2-1 e),   we  have,    again     to 
0(6), 


-1-  ♦  Vx(nJi    -  VqxB^    -  v^xBq  -  VixB^)   -  0   .  (3-2) 


3t 


We  note  that,  since  <v^xB^>  is  a  mean  quantity,  the  term  V''<Vi'<B^>  is 

smaller  than  the  term  Vx(v^xBp   and  was  dropped  in  Eq.    (3-2). 

Subtracting  Eq.    (2-6c)   from  Eq.   (2-1  c),  and  Eq.  (2-12c)  from  Eq. 
(2-1d),  yields 


9Pl 


^^  *  ^'(Pq^I  *   Pllil^  "  °  ^3-3) 


dp, 

— 1-  +  v^-VPq  +  YPqV  -v^  «  (Y-1)V.(<7Tt)  (3-^) 


to  the  lowest  order.  In  Eq.  (3"^),  we  may  write  <V  -  <,|^||  *  '^i^i' 
assuming  <  =0(1)  and  <  -  0(6^)  as  before.  Here  V  and  V  denote 
derivatives  parallel  and  perpendicular  to  the  total  magnetic  field  line 
B  +  B, ,  respectively.  Since  V  is  the  same  as  the  derivative 
perpendicular  to  the  mean  magnetic  field  B^  to  the  lowest  order,  we 
actually  define  7  as  7  -  V  -  MB^-V),  where  b  -  Bq/|Bq|^. 

By  rewriting  the  right-hand  side  of  Eq.   (3-1)  using  the  relations 
J„  -  VxB-,  and  J.  -  VxB.,  we  have 


po^^  *  yi|-'hi  -  -'(Pi*!o*Bi)  *  (§o'^^§i  *  (5r'^)5o  *  (§i|-^^li  • 


The  perpendicular  component  of  this  equation  is 
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'.(Pl    *  !o'?l^   ■  ° 


1 


to  the  leading  order.  By  assuming  the  mode  to  vanish  away  from  the 
rational   surface  as   it  does   in   the  linear  theory,    we  conclude   to  0(6) 

Pi    *  §o-§l    -  0    • 

Therefore,  the  parallel  component  of  B^  is  proportional  to  -p^ .  From 
the  equation  V«B  -  0  together  with  the  ordering  assumptions  and  the 
relation  7'Bq  -  0,  we  have  to  0(1) 


V  -B,  •  0  . 
1'^ 


2 

It  follows  that  there  exists  a  scalar  function  A  of  order  6  such  that 


B^  -  V  Axb  -  p^b  (3-5) 


Similarly,      from     Eq.      (3-'4),    we  conclude  V   'v^    -  0   to  0(1).     Therefore 

2 
there  exists  a  scalar  function  ((i   of  order  6     such  that 


v^    -  7  (J>xb  -  vb    ,  (3-6) 


where  we  write  -vb  as  the  parallel  component  of  the  velocity  v^ . 

With  these  new  functions,   the  perpendicular  component  of  Eq. 
(3-2)  becomes 


3A 

3t    -"x"  r"-    T 


■:^   +  (V,A«V,(<))b  -  nA,A  +  Bq-V(|)  ,  (3-7) 
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and      the     parallel     component     of     Eq.        (3"1),   along  with   the  relation 
?o*§1    "  "Pi '    ^^  given  as 


^  Pi    -   (5o*§l|)-'^  -  ^r^(Po*?o^    -  §o^'|-^>    -  nA  Pi    -  0   ,      (3-8) 


where  d/dt  -  3/3t  +  y.  -V,  and  A  -  7^.  Operation  B^-Vx   (Eq.    (3-1)] 
yields 


Po  ^  Y  =  (|o*§1  I  )-^(Y^  -  2bxV[pQ  *  ^  b2].7pi 


2b-V[l  b2](a  A)  *  u  A  (3-9) 


and   the  parallel  component  of  Eq.      (3-1  )  becomes 


Po  HI  ^  '    ^§o^§l|)-^Pl    *   §1|-^Po  *  ^.^.^   •  ^3-10) 


'1  '        -'1       °         11 


By  using    the  relation 


!i  -  ^  -  li 

Po  "  Po       ^o 


(3-11) 


which  holds  to  the  lowest  order,  we  eliminate  p^  from  the  system  and 
derive  the  equations  for  Pi  and  T^.  We  eliminate  the  term  ^.'^i  from 
Eqs.   (3-^*)  and  (3-8)  and  obtain  an  equation  for  p^: 
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o       to 


(Y-DB^ 

V.[(<,    V, ,    *  <   V,)T,j    .  (3-12) 


^Po*B^  II    II  11      ' 


Likewise,      from     Eqs.        (3-3),    (3-^)   and   (3-11),    we  get  an  equation  for 


f  Y  —  1     ^  T 

^^1    *   ^li-^^o    - t2Yi     .V[p^   *^  Bo^    *    (VBl|)-Vvi    *    nApJ 

-L  B^   *   YPo  1  1 


^^H-1^V.[(<      V        *<    V,)Tj    .  (3-13) 

Po       YPq   *    §2  I  I     I  I  11' 


Either  Eq.      (3-12)   or  Eq.      (3-13)   may  be  replaced  by   the   following 
entropy  equation: 


-^^    ^-U*'^'°"^'*^^'||'||    *  Vl^'lJ    •  ^3-1^) 


where  we  note  that,  to  the  lowest  order, 


y     '^1   Pi 
0   '^o 


We  now  further  simplify  these  equations  assuming  that  the  plasma 
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is  confined  in  a  cylinder  and  all  the  mean  quantities  depend  only  on 
its  radius  r.   It  is  convenient  to  introduce  the  following  independent 
variables, 


■"o  "^0°   |o| 


e  =  lole  ,         T  -  — —  t  , 


(3-16) 


VPo 


Here  (r,e,z)  denote  the  usual  polar  coordinates  and  r^  denotes  a 
particular  radius  in  the  vicinity  of  which  we  consider  the  motion  of 
the  modes.  Bq  and  B_  denote  the  azimuthal  and  longitudinal  components 
of  the  mean  magnetic  field  B^,  respectively.  We  note  that  the  6- 
direction  (with  fixed  x  and  y)  is  the  direction  along  the  equilibrium 
magnetic  field  line  ^^(r)  and  not  the  azimuthal  direction  anymore. 
Furthermore,  we  define 


B  -  |Bol  . 


M  -  — ^  ,  (3-17) 

B. 


0  --  U   . 


Where  prime  denotes  d/dr  and  every  quantity  here  is  evaluated  at 


r  -  r^. 


The  multiscale  assumption  allows  us  to  take  every  fluctuating 
quantity  to  be  a  function  of  x,  y,  6  and  t  as  well  as  r^.t,    while  every 
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mean  quantity  depends  only  on  r^  and  t.   From  the  assumption  of  the 
ordering  we  made,  we  have 


9x    3y     6 

^  -  0(1)  ,      -L  -  0(1)  ,  (3-18) 


36 


^  -  0(6^)  and  _L  -  0(1) 


We  thus  find  that  Eqs.  (3-7)-(3-10) ,  (3-12)  and  (3-13)  are  the 
nonlinear  partial  differential  equations  with  independent  variables  x, 
y,  6  and  t,  and  with  constant  coefficients  depending  only  on  parameters 
r^  and  t. 

Along  with  definition  (3-17),  we  also  introduce  the  following 
quantities 


D 


■2rp„  2p 


0 


b2o2  '  b2  ' 


4b2            rB„ 
S  -  L  .      R  -  L 


B^o^  '  n  /Pq 


.  M  - 

Y-1    d2„2   0 


Y    "2rp    ,  ^1 

—  •  ^  T„  ,      M-  1 (3-19) 


B^o'^  rB.  /r 


6  '^o 


^     Y-1   IT'S 


K.  -Ill  i 


'   rB2  /p^  '       1    ^   rBg  /T^ 
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where  every  quantity  is  evaluated  at  r  -  Tq  as  before.  The  linear 
stability  criterion  for  the  ideal  modes  is  given  by  the  inequality 
D  <  1/4,  which  is  known  as  Suydam's  criterion.  ^  We  assume  that  this 
stability  condition  is  always  satisfied  so  that  the  plasma  is  ideally 
stable  and  the  instability  arises  from  the  finite  diffusion 
coefficients. 

The  dependent  variables  are  defined  in  a  dimensionless  form  as 


1    r  0  ..  -    "^Po   r  0 


T  «  •  T,  ,    3  "  —  s, 

^-1    b2|o|3/2   1  b2|o|3/2   ^ 


where  we  have  the  relation 


(3-21) 


from  Eq.   (3-15).   Here  we  note  that  A  and  <p   are  of  order  6  and  p,  v, 
T  and  s  are  of  order  6. 

We  now  rewrite  Eqs.   (3-7)-(3-10) ,  (3-12)  and  (3-13)  in  terms  of 
the  new  variables  and  parameters  defined  above: 


£^-ii*lZA  (R-i) 

dT    _   f  1 
36 


■2H- 


_  A   4.    -  _  A    A    *    Ia.A   a|    -  IE   +   MA^^  (R-2) 

<iT  1       .  1  1       9y       1 

86 


d    -  Y6      r„   St   ^    3v         .-    -,         1    r   "i 


-   d11  .  i^llll      (K,,    Z,,    *   K   Z,]T  (R-3) 


3y  2.Y6        ^    II       II  II' 


d    -         Y6     fe  90  ^  3v       jr  -,        1    7  - 1 

36 


'f^^'^^'ii  ^1  *  v/  ^'"'^ 


dv       3d  I  ~   ~  1            3A          -   - 

—   -   ^  *  A.p      +    D   ^   +   MA   V 

di  3y           J_ 

36 


(R-5) 


Either  Eqs.      (R-3)   or   {R-^}   may  be  replaced  by    the  entropy  equation 


^"s-    (D-e)|i*    [K||    Zji    *   Y^JT    .  (R-6) 


Here  we  have  used  the  following  definitions; 
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|r,s|-<Y.Y>--i|f-|i|f 


If   .i!l*  [A^  Ha.  ill  *±  lA.f}]  *  U.  lA.f) 

II        -p      l°l 


(3-22) 


1    3x     3y'    3y2 


The  boundary  conditions  of  this  set  of  equations  will  be  discussed  in 
the  next  section. 

In  the  reduced  set  of  equations  (R),  there  are  eight  parameters, 
i.e.,  D,  0,  S,  6,  R,  M,  K  and  K  ,  which  determine  the  solutions. 
However,  if  we  apply  the  following  scale  transformations,  the 
coefficients  1/R,  M  and  K  in  those  equations  will  be  formally  replaced 
by  1 ,  RM  and  RK  ,  respectively: 


X  >  R"''/^  X  ,      y  *  R"^'^^  y 


A  +  R"""!  ,  *  -   R'\ 


p  -  R"''^^  p  ,       T  *  R"^^^  T 


V  -  R"''^^  ~         3  ^  R-1/2  - 


(3-23) 


Therefore,  using  the  following  definition  of  the  new  parameters: 
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Mr  -  MR  -  u  /nPo 

(3-24) 
Kfj  -  YK  R  -  (Y-i)<  /nPo  , 


the  solutions  of  the  system  transformed  from  Eqs.  (R)  by  the 
transformation  (3-23)  depend  only  on  Mp  and  K^  as  well  as  D,  6,  S,  B, 
and  K  but  not  on  the  magnetic  Reynolds  number  R,  so  long  as  the 
boundary  conditions  do  not  depend  on  R  after  this  scale  transformation. 
This  is  certainly  the  case  if  we  consider  the  solutions  that  decay 
rapidly  as  |x|  ,  |y|  -  ». 

The  set  of  equations  (R)  is  suitable  to  describe  the  resistive 
g-mode  fluctuations  in  a  screw  pinch  or  an  RFP,  in  which  the  azimuthal 
(or  poloidal)  field  Bq  and  the  longitudinal  (or  toroidal)  field  B  are 
of  the  same  order  of  magnitude,  and  the  plasma  beta  6,  which  is  defined 
as  the  ratio  of  the  thermal  pressure  p^  to  the  magnetic  pressure  B^/2, 
is  of  order  1.  However,  if  we  consider  a  straight  cylindrical  plasma 
with  low  6  and  a  strong  longitudinal  field  B^,  we  may  further  simplify 
the  set  of  equations  (R).  This  further  simplified  system  is  a 
reasonable  model  of  a  tokamak  plasma  with  a  large  aspect  ratio. 

Let  Eg  be  the  inverse  aspect  ratio,  i.e.,  e^  -  a/L,  where  a  and  L 

denote  the  radius  and  the  length  of  the  cylinder,   respectively.   We 

assume  that  e,  is  a  small  number  and  that  8  -  B^/B,  -  e  and  o  -  /I" 

o  z    a  a 

namely,  low  6,  strong  longitudinal  field  and  low  shear.  The  diffusion 
coefficients  such  as  the  resistivity  n,  the  perpendicular  viscosity  v 
and  the  perpendicular  heat  conductivity  k  are  also  taken  to  be  of 
order  e^.6  .  Under  these  assumptions,  it  follows  that 
D  -  0  -  R  -  M  -  K  are  of  order  1  and  S  is  of  order  e.   It  also  follows 
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that  the  safety  factor  q(r),  defined  by  q(r)  -  rB^ >) /LBg(r) ,  is  of 
order  1 ,  which  is  the  case  in  tokamaks.  By  rewriting  y  as  y 
-  |o| ^^^(o/|o| )[z/L  -  q(r)e}/q,  where  o  -  r^q'/q,  and  z/L  corresponds 
to  the  toroidal  angle,  quantities  of  order  e~  such  as  L  and  p  do  not 
appear  explicitly  in  the  equations. 

In  order  to  focus  on  the  effect  of  the  temperature  gradient,  we 
further  assume  that  the  density  gradient  is  small,  i.e.,  f^o^o^^o 
«  O(e^).  This  assumption  makes  Eqs.  (R-3)  and  (R-4)  identical  to  each 
other  to  the  lowest  order.  In  fact,  the  parameters  D  and  0  defined  in 
Eqs.   (3-19)  are  now  related  to  each  other  as 

D  -  III  e  *  i  N  ,  (3-25) 

where 


■2^oPo 


Po°^ 


is  the  parameter  related  to  the  mean  density  gradient.  Therefore,  it 
follows  that  D  -  CY-Oe/Y.  Consequently,  the  perturbed  temperature  T 
and  the  perturbed  entropy  s  are  expressed  in  p  as 

T  -  Yp/(Y-1) 

(3-26) 

i  -  p/(Y-1)  . 

Thus  the  subsidiary  expansion  of  Eqs.   (R-l)-(R-5)  in  e^  leads  to 
the  following  set  of  equations: 
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dA  a*  1  -  - 


_  A  $  -  _  A  A  +  U.A  A   -  |£  +  MAf(t)  (T-2) 

59 


A  p  .  -D  |1  *  ,l^p    .  ,  (T-3) 


Here  x  =  '^^.  and  the  equation  for  v  is  given  by  (R-5),  which  is, 
however,  decoupled  from  the  set  of  equations  (T) .  Although  we  ignored 
the  parallel  heat  conductivity  in  the  system  (T),  we  can  introduce  it 
simply  by  replacing  xA  by  X,.A  ♦  xA  in  Eq.  (T-3),  where  x..  = 
YK  .  The  boundary  conditions  to  be  imposed  on  this  system  will  be 
discussed  in  Sections  IV  and  V.  We  note  that,  as  in  the  case  of  the 
system  (R),  the  scale  transformations  (3-23)  replace  the  coefficients 
1/R,  M  and  x  in  Eqs.  (T)  by  1  ,  Mp  and  Kp,  respectively,  so  that  the 
solutins  of  the  transformed  system  depend  only  on  the  parameters  D,  Mp 
and  Kp,  if  we  impose  proper  boundary  conditions. 

These  two  reduced  sets  of  equations,  i.e.,  Eqs.  (R)  and  Eqs. 
(T),  form  the  bases  of  our  models  of  the  nonlinear  resistive  pressure 
driven  mode.  From  the  next  section  on,  we  will  investigate  these 
reduced  equations  in  order  to  characterize  the  anomalous  transport 
caused  by  these  fluctuations. 
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IV.   COUPLING  OF  THE  DYNAMO  EFFECT  WITH  THE  ANOMALOUS  HEAT 
TRANSPORT  IN  RFP  PLASMAS 


In  this  section,  we  will  consider  the  initial-boundary  value 
problem  of  the  system  of  Eqs.  (R)  and  derive  a  relationship  between 
the  dynamo  effect,  or  the  anomalous  electric  field  e  defined  in  Eq. 
(2-9),  and  the  anomalous  heat  transport  <SiV^p>:  in  fact,  we  will  show 
that  the  relationship  given  by  Eq.  (2-13)  holds  locally,  not  as  the 
integral  over  the  total  plasma  volume,  if  we  consider  the  fluctuations 
caused  by  the  resistive  g-mode  fluctuations.  Before  deriving  this 
relationship,  however,  we  will  first  discuss  the  boundary  conditions  of 
the  system  of  Eqs.  (R)  and  subsequently,  the  averaging  operations 
applying  to  the  solutions  of  this  system. 

In  order  to  make  it  easy  to  compare  the  solutions  of  the  system  to 

8  1 4 

the  more  familiar  linear  solutions,   we  transform    the  independent 

variables,  x,  y  and  e  to  the  new  variables  x  •  x,  y  -  y  +  ex  and  z  -  9. 

This  transformation  transforms  the  derivatives  in  Eqs.   (R)  as  follows: 

3/9x  »  S/3x  +  za/3y,  9/9y  -  3/ay,  3/39  -  3/3z  +  x3/3y  and  A  -  3^/3^x 

+  3  /3^y.  Although  we  will  generally  consider  the  solutions  which 

depend  on  z  as  well  as  on  x  and  y,  the  well  known  linear  solutions  of 

8        ~ 
the  resistive  fast  interchange  mode  is  the  z-independent  solution  of 

the  linearized  version  of  the  system  of  equations  (R). 

Throughout  this  thesis,  we  will  use  two  coordinate  systems, 

(x,y,9)  and  (x,y,z),  using  whichever  system  makes  equations  simpler. 

However,  it  is  easier  to  write  the  boundary  conditions  of  the  system  in 

terms  of  (x,y,z)  rather  than  (x,y,9)  because  of  the  periodicity  of  the 

z-coordinate.   From  the  definition,  the  z-direction  with  fixed  x  and  y 

is  the  direction  of  the  mean  magnetic  field  B^  at  r  -  r^.   Since  the 
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surface  r  -  r^  is  taken  to  be  a  rational  surface,  every  field  line  on 
this  surface  is  closed.  Therefore  the  particular  mean  field  line  we 
consider  on  this  surface  comes  back  to  its  original  position  after  a 
finite  length  z  -  L^.  Thus,  provided  that  this  particular  mean  field 
line  never  gets  too  close  to  itself  on  this  surface  before  coming  back 
to  its  original  position,  it  is  natural  to  impose  the  periodic  boundary 
condition  in  z  with  period  L^,  together  with  proper  boundary  conditions 
in  X  and  y.  On  the  other  hand,  the  9-direction  with  fixed  x  and  y  is 
the  direction  of  the  mean  field  line  starting  from  a  point  off  the 
rational  magnetic  surface  (if  x  *  0)  so  that  there  is  no  periodicity  in 
e  generally. 

Thus  we  impose  the  following  boundary  conditions  which  are  similar 
to  those  imposed  on  the  linearized  problem  of  the  system  (R): 


(I)     <P,    Ad),  A,  p,  V,  T  -  0      at  |x|  -  6,^ 


(II)     0,  Ad),  A,  p,  V,  T  -  0     at  |y|  -  6y 


or,  instead  of  (II) , 


(II')     (ti,   ^  <p ,   A,  p,  V,  and  T  and  their  y-derivatives 
are  periodic  with  period  26^   in  y 


and 


(III)     if),   A  0,  A,  p,  V,  and  T  are  periodic  with  period  L_  in  the 
1  ^ 


z-direction. 
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Here  6^.,  6  and  L  are  the  prescribed  positive  values  such  that 

6  <<  6vi<5,,  <<  a  and  L,/a  >  0(1),  where  6  is  the  small  scale  parameter 
X   y  i.       — 

introduced  in  Section  III  and  a  is  the  plasma  minor  radius.  6,^  and  6 
give  the  size  of  the  boundary  layer  whereas  6  give  the  size  of  the 
resistive  layer.  In  the  linear  theory  in  Ref.  8,  6,^/5  is  taken  to  be 
»,  in  which  case  the  asymptotic  behavior  of  the  linear  solution  at 
large  |x|  (i.e.,  |x|  >>  6)  takes  an  exponentially  decaying  form.  We 
presume  from  this  linear  asymptotic  analysis  that  the  solutions  of  the 
nonlinear  system  (R)  are  essentially  independent  of  the  choice  of  6j^ 
and  decay  rapidly  as  |x|  becomes  large,  as  long  as  6jj/6  is  taken  to  be 
large  enough.  Appendix  B  contains  further  discussions  of  the  boundary 
conditions. 

It  is  natural  to  define  averaging  operations  in  the  following  way 
for  a  function  f (x,y,z,T ;rQ, t)  satisfying  boundary  conditions  (I),  (II) 
(or  (II')]  and  (III): 

f  ^  dy  f  ^  dz  f(x,y,z,T;r-.t)         (U-1  ) 


<f> 


^(5„6„L,  J 


«x     -  r'y        ''- 

dx 


-6.    ■'-6, 


X  y  z   "x     "y 


f  «  lim  — 
T    T  J 


T     

dTf(x,y,z,T;rQ,t)  .  (4-2) 

0  ° 


It  is  clear  from  the  boundary  conditions  that  the  <  >  average  of  any 
spatial  derivative  of  f  with  respect  to  x,  y  or  z  vanishes.  We  also 
define  the  averaging  "<  >  as  the  combination  of  these  two  averaging, 
i.e.,  <f>  5  <f>.  It  is  easy  to  check  that  <  >,  "  and  "  >  satisfy  the 
Reynolds  conditions  (2-5),  where  the  condition  (2-5d)  holds  only  in 
terms  of  the  slowly  varying  variables  r^   and  t. 

We  point  out  that  if  we  take  ^^   to  be  large  enough  and  we  obtain 
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the  solution  f  of  the  system  (R)  that  decays  rapidly  as  |x|  becomes 
large,  then  the  average  <f>  of  the  solution  f  given  by  Eq.  (4-1)  is 
not  well  defined.  In  fact,  provided  that  such  solutions  f  depend  only 
weakly  on  the  choice  of  6^^  for  large  enough  6^^,  the  integral  in  x  used 
in  Eq.   (4-1) 


f  dx  -►     f  dx  < 

J -6^       J- 


as  6^  -»  »,  and  so  the  value  <f>  scales  as  1/6^  for  a  large  6,^.  In 
order  to  keep  the  physical  meaning  of  the  average  of  this  kind  of 
solution,  we  redefine  the  operator  <  >  as 


<f>  -  ! dx      dy 


dz  f(x,y,z;rj^,t)  .         (4-3) 
0  ° 


If  we  choose  the  boundary  condition  (II),  instead  of  (II'),  then  we 
also  need  to  modify  the  averaging  in  the  y-direction  in  a  similar  way. 
Here  we  note  the  <  >  averaging  given  by  Eq.  (4-3)  does  not  satisfy  the 
condition  (2-5)  of  the  Reynolds  conditions,  and  therefore  it  is  not  an 
averge  in  a  rigorous  sense.  However,  this  modified  average  <  > 
provides  the  contribution  of  the  fluctuating  quantity  per  unit  volume. 

In  terms  of  the  scaled  variables  defined  in  Eqs.  (3-16)  and 
(3-20)  and  the  averaging  operator  <  >  defined  above  (either  Eq.  (4-1) 
or  Eq.  (4-3)],  the  anomalous  heat  flux  "Po^^^i  ^i  r'' '  "here  the  radial 
velocity  v^^,  -  vVr,  is  expressed  as 
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B^B^o^ 


■Po<3iVir> 


3y 


(4-i4) 


On  the  other  hand,  the  averaged  electric  field  <Yi«b^>  can  be  written 
as 


Pi        V 
<v,xbi>  «  <(V  ())>«V  A)«b>b  +  <--  V  d>  -  — .  7  A> 


^iri^r 


1  1 


B^ 


1    b2  1 


/T-   3x  3y   3y  3x  -    2/p~     ^^     ^^ 


Z   M>  (VrxB_)  , 
3y      -'-' 


to  the  lowest  order.   In  particular,  for  e  -  <v^'<B^>  -  —  '^PiYi'^'lo'  ^^ 


follows 


£•§0  -  -=r   <U.*}> 


and 


£•^.0 


Po^e°   2*Y6 


2/f; 


Y6 


[-  <-p  ii>  *  ^^^-^^  <i  li>  *  j:l  <;  |^>]  .   (14-5) 

^     3y     2+Y6     3y    2+Y6    3y 


Here  we  have  used  the  relations 


■3U- 


Jo-!o.     Bo'^Po 


^0  -  i^K  ■  -T- 


3  -  _L  ^  -  _L  ^ 


We  also  note  that  U.iji}  has  a  divergence  form,  i.e., 


dx     3y    9y    3x 


so  <iA,())}>  vanishes  because  of  the  boundary  conditions.  Therefore  e'E^ 
is  of  order  6^  while  e'J^  is  of  order  6^  for  the  resistive  fast 
interchange  mode. 

We  now  confirm  Eq.  (2-11),  which  relates  the  dynamo  term  e  and 
the  anomalous  heat  flux  <SiV^^>,  directly  from  the  set  of  the  reduced 
equations  (R).  Adding  up  A  A  x  Eq.  (R-1)  and  0  »  Eq.  (R-2)  and 
taking  the  <  >  average  of  the  resulting  equation  yields  an  equation  for 
the  energy  of  the  perpendicular  magnetic  field  and  the  perpendicular 
flow 


A  ^  (<|  Y|2  *  I Yi>]  ■  -  <p  l7>  - 1  <i  Y'"' " "  'I  V''' '  ^'''^^ 


where 


7,  -  Vr(i-  -  e  ±]  ■    -°    ^ 


I  ^3x     9y'    |Bq|  9y 
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Likewise,  by  adding  up  Eqs.   (R-3),  (R-5)  and  (R-6)  with  the  following 
factors: 


p.(H-3)  .^;.(R-5)  *£^s.(R-6) 


and  by  taking  the  <  >  average,  we  obtain 


2  9t  ^  ^     2+Y6         2+Y6 


Y6 


-[D  -  _L£-  S]  <p  ^> 


2+YB 


il>  -  iillll  (e-D)  <i  ^> 

3y     2*YB         9y 


Y6  r..~    3A^    Y6   1 


2s  _   ^S 


:,i2. 


2+Y6     3y    2*Y6  R   '  J_  '      2  +  Y6    '  J_  ' 


i<I^(K||<|H-(^)|;,T|p>.Y|5^T|^>) 
36 


(4-7) 


For  the  solutions  such  as  saturated  modes  or  stationary 
turbulence,  we  can  drop  the  left-hand-sides  of  Eqs.  (U-6)  and  (4-7)  by 
taking  the  long  time  average  in  t.  Taking  the  -  average  of  Eqs.  (4-6) 
and  (4-7)  and  rearranging  the  terms,  we  obtain 


<p  il>  -  -  1  <|A  A|2>  -  M<|Ai|^> 
9y      R   '  J_  '       '  J_  ' 


(4-8) 


and 


■36- 


-D[<;  il>  - 112:11  -lh>  -Jl.  -z  ^>]  -  ^(^-^)  e  <;  i*> 

3y     2+Y6     3y    2+Y6    dy  2+Y6      3y 


2*Y6     ay    2  +  Y6  R   '  J_  '     2+Y6     '  J_  ' 


2(Y-1) 


3T 

de 


*^^^'ll"--  ^T^^  ^'-'T'  *  Y'V'''^  *  ^'"^^ 


It  is  clear  that  the  fourth  term  of  the  left-hand  side  of  Eq.  (4-9)  i3 
proportional  to  <s^v^j,>  and  the  sum  of  the  first  three  terms  is 
proportional  to  e-J^  from  Eqs.  (4-4)  and  (4-5).  By  substituting  Eq. 
(4-8)  into  <p  3(})/3y>  in  the  right-hand-side  of  Eq.  (4-9),  we  see  that 
the  right-hand-side  is  a  non-negative  function  related  to  the 
dissipation  of  the  system.  In  fact,  multiplying  Eq.  (4-9)  by 
-(B2BQo'^)(2+Y6)/4Y6rQ/p^  yields  the  following  relation: 


1  0   1 1   M       11 


Here  we  have  used 


|JJ2  .!!i:  (slyl^  .  lyl^) 


and 


'-^''■^^''Y''*  'V''^  • 
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which  are  derived  from  the   following  relations; 


1     ■         I 


We  have  thus  shown  that  Eq.   (2-13)  holds  locally,  in  other  words,   Eq. 
(2-11)  holds  on  each  rational  surface  wi th  f  -  0. 

To  conclude,  we  have  derived  in  Eq.  (4-10)  a  relation  between  the 
dynamo  effect  and  the  anomalous  heat  transport  on  each  rational  surface 
of  an  RFP  plasma,  based  on  the  system  (R)  of  the  resistive  g-mode 
fluctuations.  This  strong  coupling  of  the  dynamo  effect  with  the 
anomalous  heat  transport  explains  the  large  heat  loss  or  the  short 
energy  confinement  time  of  RFP  plasmas,  where  it  is  believed  that  the 
dynamo  plays  an  important  role  in  sustaining  the  equilibrium 
configuration  through  relaxation  processes. 
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V.   SOME  PROPERTIES  OF  THE  ANOMALOUS  HEAT  TRANSPORT  OF 
TOKAMAK  PLASMAS 


We  now  turn  to  a  model  of  the  resistive  g-mode  fluctuations  in  a 
tokamak  described  by  the  system  of  equations  (T).  In  this  section,  we 
will  show  that  the  dynamo  effect  decouples  from  the  anomalous  heat 
transport  in  tokamak  plasmas  and,  therefore,  the  anomalous  heat  flux 
transports  heat  in  the  direction  opposite  to  the  mean  temperature 
gradient. 

As  in  the  RFP  case,  the  initial-boundary  value  problem  of  this 
sytem  is  posed  more  easily  by  introducing  the  transformation  of  the 
independent  variables  x  -  x,  y  -  y  +  ex  and  z  -  9.  The  boundary 
conditions  for  A,  (f ,  A41  and  p  in  terms  of  the  new  independent  variables 
x,y,z  are  also  given  by  (I),  (II)  (or  (II')]  and  (III)  in  the  previous 
section.  (See  also  Appendix  B.)  The  averaging  operators  <  >  and 
defined  in  the  previous  section  are  also  used  in  this  section. 

The  relation  between  the  dynamo  dissipation  g'Jq  and  the  anomalous 
heat  transport  "Pq'^SiX^  given  by  Eq.  (4-10)  for  an  RFP  plasma  is 
considerably  simplified  for  a  tokamak  plasma.  Since  s  -  p/(Y-1 )  [Eq. 
(3-26)]  for  the  tokamak  case,  Eq.   (4-5)  becomes 


to  the  lowest  order  under  the  expansion  in  the  inverse  aspect  ratio  e^. 
From  the  coefficients  of  the  right-hand-side  of  Eq.  (5-1),  it  is  clear 
that  e*J-  is  of  order  e^.  Here  we  recall  that  we  measure  the  length  in 
r^  units  and  the  magnetic  field  in  B  units,  so  that  from  the  assumption 
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of   the   Eg-expansion   discussed   in   Section   III,    we    have 

p  -  Bg  -  0  -  Eg.   On  the  other  hand,   the  radial  component  of  the 

2 
anomalous  heat  flux  "Po'^^i^ir^  given  by  Eq.    (U-4)   is  of  order  e^, 

which  is  of  the  same  order  as  the  collisional  heat  transport  <|Tq. 

Thus  the  heat  production  e'Jq  related  to  the  dynamo  effect  is  smaller 

by  an  order  of  the  inverse  aspect  ratio  e,  than  the  anomalous  heat 

transport. 

Taking  into  account  the  ordering  in  e^     and  recalling   the 

assumption  that  K   does  not  enter  the  equations,  we  can  simplify  Eq. 

(U-9)  as 


-D<;  |1>  -  x<|vi|2>  .  (5-2) 


Here  Eq.  (3-26),  and  the  relations  8  -  YD/(Y-1  )  and  YK  -  x  are  used: 
we  note  that  the  anomalous  heat  flux  -Po^^i^ir''  ^^  proportional  to 
<p3(ti/8y>  in  the  tokamak  case.  We  may  also  obtain  Eq.  (5-2)  in  the 
following  way:  by  multiplying  Eq.  (T-3)  by  p  and  taking  its  <  > 
average,  we  get 


li.<|;|2>   .   -D<p|l>   -   x<lYl'>    •  (5-3) 


Taking  the  ~  average  of  this  equation  yields  Eq.   (5-2). 

In  terms  of  the  physical  variables,  Eq.   (5-2)  may  be  written  as 
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1 


which  corresponds  to  Eq.  (4-10),  where  we  recall  our  assumption  that 
Pq  "  Pq'^o  ^°^  the  system  (T).  Here  we  used  Eqs.  (3-16),  (3-20)  and 
(3-26). 

To  conclude,  we  have  shown  that  the  anomalous  heat  transport  is 
decoupled  from  the  dynamo  effect  in  tokamak  plasmas.  We  have  also 
shown  from  Eq.  (5-4)  that  the  anomalous  heat  flux  across  the  magnetic 
surface  is  in  the  direction  opposite  to  the  mean  temperature  gradient. 
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VI.   NONLINEAR  STABILITY  ANALYSIS  OF  A  TOKAMAK  PLASMA 

In  this  section,  we  will  discuss  the  stability  of  the  null 
solution  of  Eqs.  (T) .  The  boundary  conditions  to  be  used  in  this 
section  are  (I),  (II)  (or  (II')]  and  (III)  given  in  Section  IV  and  V. 
The  null  solution  of  this  system  represents  the  equilibrium  state  of 
the  plasma.  We  define  the  stability  of  the  null  solution  in  the 
following  way:  let  ed)  be  the  energy  norm  of  the  system  defined  as 


e(T)  -  ^  <|V  A|2  *  ivil^  *  |P|^> 


Although  this  ed)  is  not  the  total  energy  of  the  perturbed  plasma,  we 
take  ed)  to  be  a  measure  of  the  amplitude  of  those  perturbed 
quantities.  The  null  solution  of  the  system  (T)  is  called 
asymptotically  stable  (in  the  sense  of  Liapunov)  if  there  is  a  positive 
number  6j  such  that  for  any  initial  condition  the  energy  norm  of  which 
is  less  than  6j,  i.e.,  for  e(0)  <  6j,  the  energy  norm  approaches  zero 
for  large  t,  i.e.,  e(T)  +  0  as  i  -»•  •.  In  particular,  if  6^  is  finite, 
it  is  called  conditionally  stable  and  if  6j;  -  »,  then  it  is  called 
globally  stable.  The  globally  stable  null  solution  is  also  called 
monotonically  stable  if  de(T)/dT  <  0  for  all  t  >  0. 

The  various  notions  of  stability  defined  above  are  based  on  the 
nonlinearity  of  the  system.  For  a  perturbation  with  infinitesimal 
amplitude,  however,  we  have  the  spectral  problem  of  the  linearized 
system  of  Eqs.  (T).  Assuming  that  the  solutions  of  the  linearized 
system  depend  exponentially  on  x  as  e^"^,  where  q  is  a  complex  number, 
we  replace  3/3t  by  q.  Provided  there  exist  such  numbers  q  for  which 
the  linearized  system  has  nontrivial  solutions,  the  numbers  q  are 
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called  the  eigenvalues  of  the  system.  The  linear  stability  of  the  null 
solution  is  then  defined  in  connection  with  the  eigenvalues  q  as 
follows:  the  null  solution  is  called  linearly  stable  if  there  are  no 
eigenvalues  such  that  Re  q  >  0;  marginally  stable  if  there  is  at  least 
one  eigenvalue  with  Re  q  -  0  and  all  the  other  eigenvalues  have 
Re  q  £  0;  and  unstable  if  at  least  one  eigenvalue  has  Re  q  >  0. 

In  Section  VII,  we  will  discuss  the  case  where  the  null  solution 
of  the  system  (T)  is  linearly  unstable.  A  small  perturbation  given  to 
such  a  system  initially  grows  with  a  growth  rate  q  and,  possibly, 
saturates  eventually  with  finite  amplitude.  In  order  to  look  for  such 
solutions  of  the  system,  it  is  useful  to  know  when  the  system  becomes 
linearly  unstable.  The  following  proposition  asserts  that  D  _<  0  is  a 
sufficient  condition  for  linear  stability.  In  other  words,  the  system 
may  be  linearly  unstable  only  when  D  >  0. 

PROPOSITION  U   Suppose  R,  M  and  x  are  all  positive  constants.   If 
D  £  0,  then  Re  q  <  0. 

PROOF:  The  linearized  equations  of  the  system  (T)  are 


qA  -  ^  +  1  A,A  (6-la) 

-   R  1 
36 


qA,^  -  -L  A,A  -  i£  +  MA^i  (6-1  b) 

1       .   1      9y       1 

36 
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qp   -    -D  |1   *    xA.p    .  (6-1c) 


Here  we  allow  A,  ip,  p  and  q  to  take  complex  values.  Adding  up 
A  A  X  Eq.  (6-1a)*  and  ip*  x  Eq.  (6-1b),  and  taking  the  <  >  average  of 
the  resulting  equation,   we  obtain 


q»<|V   A|2>    *    q<|Y|2>    -    -    <;  ^>    -  ^   <|  Y|2>    -   M<|  y|2>    .       (6-2) 


Here     *     denotes     the  complex  conjugate.     Similarly,   by  multiplying  Eq. 
(6-1 c)*   by  p  and  taking   the  <    >  average,    we  have 


q*<|p|2>   =   -D<p  ^>   -   x<|V,p|2>    .  (6-3) 


If  D  -  0,    then  the  real  part  of  Eq.      (6-3)   becomes 


(Re   q)<|p|2>   -   -x<|V   P|^>    . 

If   <|V  p|^>   »  0,    then  Re  q   <  0.      If   < | V  p| ^>   »  0,    or  p  5  0     because     of 
the  boundary  conditions,    then  the  real  part  of  Eq.      (6-2)   becomes 


(Re   q)(<|7A|2>    +    <|V,i|^>)    -   --1  <|AA|2>    _   M<|A4ir>    • 
1  1  "1  1 


Therefore  Re  q  is  negative  for  any  nontrivial  solutions  of  the  system 
(6-1).  If  D  <  0,  then  eliminating  the  term  <pa(t>«/9y>  from  Eqs.  (6-2) 
and    (6-3)   and   taking   the  real  part  of   the  resulting  equation  yields 
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-(Re  q)(<|V  A|2>  *  <|Y|2>  -  ^  <|;|2>] 


Therefore,  for  the  nontrivial  solutions  of  the  system  (6-1)-(6-3),  we 

have  Re  q  <  0. 

(Q.E.D.) 

We  now  turn  to  the  nonlinear  modes.  In  previous  sections,  we 
discussed  some  properties  of  the  solutions  that  satisfy  the  conditions 
that  the  "  averages  of  3<|7  A|2>/3t,  9<|V^|^>/3t  and  3<|p|^>/3t 
vanish.  The  saturated  modes  and  stationary  turbulence,  which  we  are 
mainly  interested  in,  are  the  "nontrivial"  special  cases  of  solutions 
with  such  a  property:  <|V  A|2>  <|V  d)|^>  and  <|p|^>  are  bounded  in  time 
T  (so  ed)  is  bounded  in  time  x  and  deTrTTdT  -  0),  but  none  of  them 
vanish  at  t  ■»•  +".  However  this  kind  of  solutions  may  not  exist  for  a 
certain  range  of  parameters  of  R,  M,  x  and  D.  In  fact,  the  folx^wing 
proposition  asserts  that  if  the  heat  conductivity  x  i3  taken  to  be  0, 
then  there  is  no  such  saturated  mode  nor  stationary  turbulence.  This 
shows  the  importance  of  the  heat  conductivity  x  to  sustain  the 
anomalous  heat  transport  for  a  significantly   long  time  scale. 

PROPOSITION  2:  Suppose  R  and  M  are  positive  constants,  x  "  0  and  D  »  0. 
If  ed)  is  bounded  in  time  t,  then  <|V  A|^>,  <|V  i|^>  and  <p34i/3y>  ■►  0 
as  T   -►  •• 
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PROOF:  Before  proving  this  proposition,  we  recall  that  <p9(|i/9y>  is 
proportional  to  the  anomalous  heat  flux.  Therefore  this  proposition 
asserts  that  the  ~  average  of  the  anomalous  heat  flux  is  zero  if  x  "  0« 
From  now  on,  we  again  consider  real-valued  solutions  A,  41  and  p  of  the 
system  (T).  From  the  assumption,  the  ~  averages  of  9<|7  A|  >/3t, 
3<|7  $|^>/3t  and  3<|p|^>/3t  vanish,  and,  therefore,  we  have  Eqs.  (4-8) 
and  (5-2).   If  x  -  0  then  we  obtain  <p  ii>  -  ~\l J\^>   -  <|a  *|^>  -  0. 

^y  1  1 

Since  <|a  A|^>  and  <|a  4i|^>  are  positive  functions  of  t,  therefore, 
<|aa|^>,    <|a^|^>   -  0   as  T   -►  0. 

In  order  to  prove  that  <|V  Ap>  and  <|V  4i|  >  -►  0,  we  will  use  the 
Poincare's  inequality :^^  for  any  real-valued  function  f  (with  a  proper 
differentiability  condition)  satisfying  the  boundary  conditions  (I), 
(II)    (or   (II' )j   and   (III),    there  exists  a  positive  constant  C  such  that 

<|V  f|^>   <   C<|a  f|^>    . 

Here  we  note  that  the  constant  C  is  finite  since  the  domain  on  which  f 
is  defined  is  finite.  Applying  the  Poincare's  inequality  to  A  and  0 
and  using  the  fact  that  <|A  A|^>  and  <|a  i|^>  +  0  as  t  -<•  »,  we  obtain 
<|V  A|^>  and  <JV  i|^>  -►  0  as  t  ■►  +».  To  prove  that  <p3<{i/3y>  -►0  as 
T   -►  0 ,    we  use   Schwartz's   inequality: 


<p|i>i^<<i;i^><i||i'> 


-ne- 

Since  <|p|  >  is  assumed  to  be  bounded  in  time  t,  <p  3()i/3y>  *  0  as 
T  ->•  0. 

(Q.E.D.) 

From  Proposition  1  and  Proposition  2,  we  find  that  we  need 
positive  R,  M,  x  and  D  in  order  to  look  for  the  solutions  A,  (>  and  p 
that  grow  from  small  initial  values  and  eventually  saturate  with  finite 
amplitude  or  lead  to  stationary  turbulence.  In  the  following  theorem, 
we  will  show  that,  if  the  plasma  is  sufficiently  viscous  or  heat 
conductive,  the  null  solution  is  still  stable  even  for  some  positive  D. 

THEOREM  j_:  Suppose  R,  M  and  x  are  all  positive  constants.  If  either  M 
or  X  is  large  enough,  there  is  a  positive  critical  value  D^  of  D  such 
that  for  any  D  <  D^,  the  null  solution  is  monotonically  stable. 

PROOF; 

Let  H  be  the  set  of  all  real  valued  functions  (with  sufficient 
smoothness;  for  a  more  specific  definition  of  H,  see  Appendix  C)  of  x, 
y  and  z  satisfying  the  boundary  conditions  (I),  (II)  [or  (II')]  and 
(III)  given  in  Section  IV.  The  solutions  A,  0  and  p  of  the  system  (T) 
are  obviously  some  elements  of  H.  For  any  element  f  in  H,  we  have  the 
following  Poincare's  inequalities  as  before:  there  exists  a  constant  C 
such  that 


<f^>  <  C<|7  f|2> 

1  (6-5) 

<|V  f|2>  <  c<|A  f|^>  . 
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By  adding  Eqs.      ('4-6)   and    (5-3),    we  obtain 


-^  e(T)    -  -(UD)    <p  |1>  -  F(t)    , 
dT  3y        - 


(6-6) 


where 


:i2s 


F(t)    -  1  <|A   A|'^>    *   M<|A.().|^>    *   X<|V,Pr> 
*^  1  1  1 


Therefore 


<-p|l> 
■^  e(T)    -   -F(t)    [l    *    (1*D)      ^/^    ] 


dT 


F(t) 


Let  us  define  D^  as 


sup 
A,0,p  e  H 


<p|l> 


F(t)    '    -   1    +    D, 


(6-7) 


~  9* 
Since  | <p  — >/F|  is  bounded  from  above,  we  have  D^  >  -1 .   In  fact,  this 

boundedness  is  shown  as  follows.   We  have 


<"^  l^>  ^ 


F(t) 


<-p  |1>2 
3y 


{M<\1  l\^>  *   x<|v  p|^>)^ 


(6-8) 


From  Schwartz's  inequality,  i.e. 


-i<8- 


<P||>'  <<|P|'>  < 


3$ 

3y 


">  <  <|p|2><|V  0|2> 


and  Poincare's  inequalities  (6-5)  for  p  and  V  *,  the  right-hand-side  of 
inequality  (6-8)  is  less  than 


<|V  p|2><|;   ^|2> 
:2 I 1 

(M<|A   0|2>    *    j^<|v  ;|2>)2 


2xM 


It  follows  that  the  supreme  value  given  by  Eq.  (6-7)  can  be  less  than 
1  by  taking  either  x  or  M  to  be  large  enough.  Therefore,  with  such  a 
choice  of  X  and  M,  we  obtain  D^,  >  0. 

Continuing  our  proof,  we  get  from  Eq.   (6-7) 


±.<„  <-,<„(,- J^j 


(6-9) 


Let  A  -  min[-,  M,x),   then,  by  using  Poincare's  inequalities  (6-5)  for 

H 

V  A,    V  A  and  p,    we  have 


F(T)    -    A(±   <|Z^A|2>    .   M   <|  Y|2>    .  I  <|  Y|2>] 

>  A(<|a   A|2>    ♦    <|A   ♦|2>    +    <|7   p|^>) 

-^  ^L  Jb 

>  AC'^-ed)    . 


(6-10) 


From  inequalities  (6-9)  and  (6-10),  we  obtain 
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^         I     \     ^  A   -  , r ,     1+D  ^ 

dT      -   C         1+D„ 


or,  by  integrating,  we  get 


e(T)  <  e(0)  exp(-^  (1  -^)tJ 

c 


Hence,  if  D  <  D^,,  then  the  perturbation  decays  exponentially  with  a 
decay  constant  proportional  to  the  smallest  dissipation  constant  A 
-  min(1/R,M,x).  (Q.E.D.) 

Since  |  <p9(ti/3y  >  /F|  is  bounded  from  above,  there  exists  the  least 
upper  bound  of  |  <p3<}i/3y>/F|  .  Suppose  A^,  (ji^  and  p^  take  the  maximum 
value  of  -<p3())/3y>/F.  (We  note  that  A^,  $q  and  p^  may  not  have  enough 
smoothness  but  they  can  be  approximated  by  functions  in  H  as  "closely" 
as  possible.  See  Appendix  C. )  By  choosing  these  A^,  <^^  and  p^  as  an 
initial  condition,  we  have,  from  Eq.   (6-7), 


^.(o,.-no,(,  -J^) 


It  follows  that  de(0)/dT  >_  0  if  D  _>  D^.  Therefore  the  following 
corollary  holds: 

COROLLARY;  Under  the  condition  of  Theorem  1,  define  D^,  by  Eq.  (6-7). 
Then,  D  <  D  is  the  necessary  and  sufficient  condition  for  monotonic 
stability. 
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We  call  this  critical  value  D^,  the  energy  stability  limit.  We 
also  define  the  linear  stability  limit  D^  as  the  value  that  makes  the 
null  solution  of  the  system  linearly  marginally  stable.  The  following 
theorem  gives  the  relationship  between  the  energy  stability  limit  D^ 
and  the  linear  stability   limit  D^. 

THEOREM  2:   Suppose   R,      M     and     x     ^'^^     ^11     positive     constants.        Then 
°c  1  Dl- 

PROOF:      Adding   up  Eqs.      (6-2)   and   (6-3)   and  taking   the  real  part  of  the 
resulting  equation,    we  obtain 


0   =    (1*D,  )    Re   <p  iil>    *  I   <|A,A|2>    +   M<|A,i|2>    +    X<|p|^>    . 
^  Sy  R  I  1 


where  we  used  the  condition   that   Re  q  »   0   if  D  =  D^^      Therefore  we  have 


1  "  3y     -^  9y 


1+D 


(6-11) 


1  <|A,A|2>  +  H<\ll\^>    *    x<|V,p1^> 


Here  we  write  p  -  Pp  +  ip^  and  $  -  4)^  +  i^^.  Writing  <|A  4i|  >  • 
<|a  ip|>  +  <|A  ij2>  and  <|V  p|^>  -  <|V  Pp|^>  *  <\V  'p^\^>  and  applying 
the  following  inequality  to  Eq.   (6-11) 
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a 
<  max    ( — ]    , 


where  C  and  bj^  are  positive,   we  obtain 


1        .                                                 9y 
-_-_   <      sup  1— 


()),p  G  H       m<|a^|2>   ♦    x<|V|P|^> 


sug 
A,(t>,p  G  H 

1 


F{t) 


1>D, 


It  follows  that  Dg  <  D^.  (Q.E.D.) 


To  conclude  this  section,  it  is  found  that  we  need  positive 
parameters  of  D  as  well  as  R,  Xf  M  in  the  system  (T)  in  order  to  find 
the  solutions  that  grow  from  small  initial  values  and  eventually 
saturate  with  finite  amplitude  or  lead  to  stationary  turbulence.  It  is 
also  found  that,  with  a  proper  choice  of  the  viscosity  M  and  the  heat 
conductivity  x.  there  exists  a  positive  energy  stability  limit  D^, 
which  is  smaller  than  or  equal  to  the  linear  stability  limit  D^.   This 
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shows  that,  if  the  plasma  is  sufficiently  viscous  or  heat  conductive, 
the  null  solution  is  stable  (to  any  perturbation)  even  for  some 
positive  D.  Since  D  is  the  only  free  energy  source  of  the  system,  if 
D  <  D^,  then  all  the  energy  fed  to  the  modes  by  the  mean  pressure 
gradient  dissipates  through  the  collisional  diffusion,  and  not  through 
convection.  When  D  >  D^^,  then  even  an  infinitesimal  perturbation  given 
to  the  null  solution  starts  to  grow  and  the  free  energy  produced  from 
the  mean  pressure  gradient  is  transfered  by  convective  motion  of  the 
plasma.  In  the  next  section,  we  will  consider  the  behavior  of  such 
convection  near  the  linear  stability  limit  D^,. 
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VII.   BIFURCATION  ANALYSIS  OF  A  TOKAMAK  PLASMA 

In  this  section,  we  will  consider  the  nonlinear  behavior  of 
perturbations  of  a  linearly  unstable  null  solution  and  we  will  derive 
the  dependence  of  the  anomalous  heat  transport  on  the  prameter  D,  v^iich 
is  proportional  to  the  mean  pressure  gradient.  As  discussed  in  the 
previous  section,  if  the  parameter  D  is  larger  than  the  linear 
stability  limit  D^,  then  modes  on  the  rational  surface  begin  to  grow 
even  for  an  infinitesimal  initial  perturbation.  It  is  possible, 
however,  that  the  given  perturbation  eventually  saturates  with  a  finite 
amplitude.  In  this  case,  the  sum  of  the  null  solution  and  this 
saturated  mode  gives  us  another  equillbrum  state,  a  stable  one,  which 
bifurcates  from  the  unstable  null  solution.  It  is  also  possible  that  a 
small  perturbation  given  at  t  -  0  developes  into  turbulence  if  the 
pressure  gradient  D  is  sufficiently  larger  than  the  linear  stability 
limit  D^.  In  the  following,  we  consider  the  former  case:  D  is  slightly 
larger  than  D^^  so  that  we  expect  to  have  steady  convection  within  the 
boundary  layer,  which  is  analogous  to  the  Benard  convection  in 
conventional  fluid  dynamics. 

For  simplicity,  we  will  assume  that  mode  rational  surfaces  in  a 
plasma  are  well  separated  from  each  other  and  mode-mode  Interactions 
between  two  different  magnetic  surfaces  are  ignorable.  This  is 
certainly  not  the  case  for  a  turbulent  plasma  but  this  is  also  a 
reasonable  assumption  for  a  certain  region  of  a  relatively  quiet 
tokamak  plasma.  Under  this  assumption,  we  will  consider  nonlinear 
modes  on  a  single  rational  surface,  where  the  modes  with  the  same 
helicity  interacts  with  each  other  nonllnearly.  These  modes  on  a 
single  rational  surface  can  be  described  as  the  z-independent  solutions 
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of  the  nonlinear  system  (T),  where  x,  y  and  z  are  the  coordinates 
defined  as  x  -  x,  y  -  y  +  ex  and  z  -  6  in  terms  of  coordinates  (x.y.e) 
used  in  this  system  (see  Sections  IV  and  V,)  We  will  use  the  coordinate 
system  (x,y,z)  throughout  this  section. 

The  z-independent  normal  mode  solutions  of  the  linearized  system 
of  Eqs.   (T)  have  the  forms 

A(x,y)  -  A(x)  exp  (imky) 

(p(x,y)  -  ip(x)   exp  (imky)  (7-1  ) 

P(x,y)  -  p(x)  exp  (imky)  , 

if  we  impose  the  boundary  condition  (II'),  instead  of  (II),  given  in 
Section  IV.  Here  m  is  an  integer,  and  k  -  Tr/6y  (6  is  defined  in 
Section  IV).  The  helicity  of  these  modes  is  given  by  the  pitch  length 
u  in  Eqs.  (3-17)  at  this  magnetic  surface.  In  fact,  we  have 
y  «  (z-ye)  to  the  lowest  order,  where  9  and  z  are  the  azimuthal  (or 
poloidal)  direction  and  the  longitudinal  (or  toroidal)  direction  of  the 
cylindrical  plasma,  respectively.  Since  6y  is  a  quantity  of  the  order 
of  1//r,  where  R  is  the  magnetic  Reynolds  number,  we  note  that  the 
modes  given  in  Eqs.   (7-1)  with  m  of  order  1  are  the  modes  with  short 

Q 

wavelengths,    which  correspond  to  the  resistive  fast   interchange  modes. 

The  z-independent  solutions  of  the  system  (T)  are,  therefore,  a 
nonlinear  extension  of  the  linear  resistive  fast  interchange  modes.  In 
particular,  the  z-independent  steady  solutions  are  the  extension  of  the 
linear  marginal  modes  of  the  system.  We  now  consider  such  steady 
solutions  on  a  single  rational  surface.  When  the  parameter  D,  which  is 
the   free  energy  source  of   the  modes,    is  slightly   larger   than  the  linear 
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stability  limit  D^,  a  perturbation  growing  from  an  infinitesimal 
initial  value  is  expected  to  saturate  with  small  amplitude.  In  this 
case,  we  are  able  to  obtain  the  expressions  of  the  nonlinearly 
saturated  modes  and  to  estimate  the  anomalous  heat  transport  caused  by 
these  modes,  by  expanding  the  solutions  in  terms  of  the  small 
amplitude, ^  ^ 

Let  t     be  the  energy  of  the  perpendicular  flow  of  such  a  solution, 
i.e. , 


e^  "  I   <l^|*l^>  •  (7-2) 


Assuming  that  e  is  a  small  number,  we  introduce  new  dependent  variables 
A,  (|)  and  p  defined  by  A  »  eA,  ip  -  t<t>  and  p  -  ep.  Since  we  consider  the 
steady  z-independent  solutions,  we  have  3/9t  -  0  and  9/3z  -  0  (i.e., 
9/9e    «  x3/3y).     Therefore,    Eqs.      (T)  become 


e{0,A}   -   X  —  <})+-  A  A  (7-3a) 

3y 


t{l,l   ^}  -  X  —  A  A  +  e{A,A  a}  -  ^  +  hI^a  (7-3b) 

3y  3y 


:i;,;}  -  -  dI*  +  xA  P  .  (7-3C) 

3y 


Our  aim  in  this  section  is  to  solve  the  system  (7-3),  together  with  Eq. 
(7-2),  assuming  that  e  is  a  small  quantity.   As  boundary  conditions,  we 
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take  (II'),  instead  of  (II),  as  well  as  (I)  and  (III)  given  in  Section 
IV,  so  that  the  solutions  of  this  system  are  the  natural  extensions  of 
the  linear  marginal  modes. 

In  the  system  (7-3),  all  the  coefficients  depend  on  e 
analytically:  in  fact,  some  coefficients  do  not  depend  on  e  and  the 
others  are  just  e's.  Therefore  the  solutions  A,  0  and  p  of  this  system 
depend  on  e  analytically  if  they  exist  at  all.  However  the  existence 
of  the  nontrivial  solutions  of  this  system  is  not  guaranteed  for  all  e 
(for  fixed  other  parameters  D,  R,  M  and  x)-  In  fact,  if  e  -  0,  then 
the  system  (7-3)  forms  an  eigenvalue  problem  and  that  D  -  D^  is  the 
condition  that  the  nontrivial  solutions  of  this  linear  system  exist. 
Thus,  for  a  finite  e,  we  choose  a  proper  D  -  D(e)  as  a  function  of  e 
(whereas  the  other  parameters  R,  M  and  x  remain  as  given  constants)  30 
that  we  expect  the  solutions  of  the  system  (7-3)  to  exist  for  any 
(small)  e.  Here  we  assume  the  analytic  dependence  of  D  on  e  and  expand 
the  solutions  and  D  in  power  series  of  e  as  follows: 


n-0  n-0  n-0 


and 


D  -  Dl  *  I     D^e" 
n-1 


We  now  seek  formal  solutions  of  this  type  by  substituting  these 
quantities  into  Eqs.  (7-3).  Considering  the  case  where  \>  ^  (this 
choice  is  always  possible  as  shown  In  Section  VII),  we  obtain  the 
following   system  as  the  zeroth  order  system  in  e: 
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where  L  is  a  linear  operator  defined  as 


(7-4) 


-2  ^  ^o       ^Po         -2-    ■ 


9y  3y 


34i 


o  _  _x 

D 


3y 


T^A   p^ 
L     1 


(7-5) 


and  Aq  is  determined  from 


Yo  -  -R^ 


3*. 


3y 


(7-6) 


The  lowest  order  term  of  Eq.   (7-2)  determines  the  amplitude  of  (p. 


<iYoi  >  ■  2 . 


(7-7) 


We  now  consider  marginally  stable  solutions  having  the  following 


forms ; 


Aq  -  Aq^ (x)  cos  ky  , 

*o  ■  ♦oi^^^  ^^"  "^y  • 

Po  "  Poi^''^  °°^  "^y  ' 


(7-8) 
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where     k     is  a  given  wave  number  in   the  y-direction.      In  this  case,   Eq. 
(7-'^)  becomes 


^01 


^01 


Rk^x^   *   MO^-k2)2 


-JLo2-k2) 


foi 


Poi 


(7-9) 


^^re  3j^  -  8/9x,  and  the  functions  (Jiqi  and  Pq^  are  taken  to  be  real- 
valued  functions  of  x.  We  think  of  L,^  as  an  operator  applied  to  a  pair 
of  real-valued  functions  (ui,U2)^  of  x,  where  t  denotes  the  transpose 
of  the  vector  and  u^  and  U2  are  sufficiently  smooth  L^-f unctions.  We 
define  the  inner  product  of  such  functions  (u^,U2)''  and  (v^,V2)^  as 


(u^ ,U2) 


l.»2J 


<u^v^  +  U2V2> 


It  is  clear  that  L,^  is  a  Hermitian  operator  with  respect  to  this  inner 
product. 

In  Eq.  (7-9),  we  note  that  the  linear  stability  limit  D^  is  a 
function  of  k.  In  order  not  to  generalize  the  problem  too  much,  we 
here  choose  the  wave  number  k  such  that  the  mode  in  Eqs.  (7-8)  is  the 
only  marginal  mode,  in  other  words,  all  the  other  modes  whose 
y-dependence  is  given  by  either  sin  X,ky  or  cos  Iky  (£,  •  1)  have 
negative  linear  growth  rate  when  D  -  D,  .  Once  we  choose  such  a  wave 
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number  k,  then  we  fix  this  k  and  D^,  and  treat  them  as  constants 
throughout  our  discussion.  In  this  case,  the  operator  L^^,  where  i  is 
an  integer,  is  always  invertible  unless  I  -  1,  so  that  the  linear 
equation  L2^,^(u^,U2)  -  (0,0)^  only  gives  the  trivial  solution  if  H  «  1. 
Thus,  by  solving  the  linear  homogeneous  equation  (7-6)  and  Eq. 
(7-9)  with  the  normalization  condition  (7-7),  we  obtain  the  solutions 


A  ,  (J)  ,  and  Pq  [in  the  form  of  Eqs.   (7-8)  j  of  the  zeroth  order  syst 


em 


(7-4).   With  these  solutions,   the  first  order  system  in  e  from  Eqs. 
(7-1  )  and  (7-3)  becomes 


<V  Oq-V  4^>   -   0 


(7-10) 


A  A 


l^cMo^    -  Rx  A   Uq,Aq1 

9y 


Uo-Mol 


(7-n) 


3y 


A  A 


A  A,    -   -Rx  —  <^^    *   RUq'^o^    • 


(7-12) 


If  we  expand  the  solutions  A^ ,  0^  and  p^  in  Fourier  series  in  y  such  as 

A  A  A  A  A  A 

A,      -  Ij^A^j^   cos   S.ky,      <p^      -  Ii*Ti   sin  iky  and  p,    -   liP-n   cos  iky,    it   is 
clear   that  we  obtain  only   three  sets  of  equations     for     £  -  0,1      and     2 
from  Eqs.      (7-11 )   and    (7-12). 
For   J,  -    1 ,    we  have 
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'11 


D,    . 


.PnJ         L"  r  "^oi^^). 


(7-13) 


The  necessary  and  sufficient  condition  that  Eq.  (7-13)  have  a 
nontrivial  solution  is  that  the  right-hand-side  is  orthogonal  to  the 
null  space  of  the  adjoint  operator  Lj^  of  L,^.  This  condition  determines 
D^ .  Since  L^^  •  L|^,  we  only  have  to  require  that  the  inner  product  of 
((foitPoi)   ^^'^  the  right-hand-side  of  Eq.      (7-13)   be  equal   to  zero: 


0    . 


(7-1U) 


Multiplying   the  second  equation  of  Eq,      (7-9)   by  Pq^    and  taking    the  <   > 
average  of   the  resulting   equation  yields 


k<Poi<Coi> 


JL  <(9,Poi)^  *  ^W\> 


so  <Poi*oi'*   ^^  nonzero  for  nonzero  Pq^  .      Therefore,    from     Eq.        (7-1^), 
D^    •  0  and  the  solution  of  Eq.      (7-13)   becomes 


'11 


'11 


^01 


POI 


(7-15) 


with  some  constant  C.   For  £,  -  0  and  2,  we  obtain  ('I'lO'PlO^  ^"^ 


^*12'Pl2^  from  Eq.   (7-11)  by  Just  inverting  the  opertors  L^  and  L2|^ 
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since     the     null     spaces  of   these  operators  consist  of  only  one  element 
(0,0)^.      By   taking   into  account  Eq.      (7-10),    i.e., 


<YorYii>  ■  0  ' 

we  have  C  -  0  in  Eq.   (7-15)  and  the  solution  of  Eq.   (7-10)  is  given 
by 


"*1' 

*10 

- 

+ 

A 

^ 

J' I 

_PlO 

<p^2   sin  2ky 


p.p  cos  2ky 


(7-16) 


The  second  order  system  in  e  from  Eqs.   (7-3)  is 


^*2i 


(7-17) 


and 


A  A2  -  -Rx  —  <P2  *  ^ 


_3_ 
3y 


where 
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f   -    Uo-Y^}    *    Ui.Yof    -    lAo-Ml^    "    ^^tVo^    -  ^  4r  h 


3y 


g    ■     Uo-Pll     *     Ui.Pq} 


and 


mm  A 

We     solve  Eq.      (7-17)   by  using    the  following   Fourier  expansion  in  y:   A2 

-  \i  k2i   cos   Iky,        <s>2       "   Ij,   *2S,   ^^"   '^'^y*        P2        "   ^^   P2i   '^°^   '^^^ '        ^ 

-  ^j^  fj_  sin  £ky  and  g  -   [^  g£   cos   iky.      For   i  -    1,   we  have 


'21 


L^2lJ 


_-(gl*kD2<J.oi)/V 


(7-18) 


the  solvability  condition  of  which  determines  D2.  similar  to  the  case 
of  Eq.  (7-13).  The  necessary  and  sufficient  condition  that  Eq. 
(7-18)  has  a  solution  is,  as  before,  that  the  inner  product  of 
(((loi'Poi^  and  the  right-hand-side  of  Eq.  (7-18)  be  equal  to  zero. 
Therefore,    D2  is  given  by 


D2   -    (Dl<0o1^1>    -    <P0l8l>^/'<<P01*01> 


(7-19) 


Assuming  that  this  D2  is  not  zero,  we  are  able  to  relate  the 

anomalous  heat  transport  <SiV^j,>  to  the  parameter  D.   Since  <p3i^/9y> 

2  «  ~  P 

-  e  k<pQ^((iQ^>  and  D  -  D^^  +  e'D2  to  the  lowest  order  of  e,  we  have 
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<p  _£>  -  i«poi0oi>  — n-^ 

z 
3y 

or  op 

BqB^o'^  <.      «  (D-D.) 

-  Pq   <3,v,^>   -       '<<Poi<|.oi>  —5 (7-20) 

(Y-1)  /Tq  ^2 

near  D  -  D^^.  Here  we  used  Eqs.  (3-20)  and  (3-26)  to  relate  <p90/3y> 
to  <s^v^^>.  Equation  (7-20)  gives  the  dependence  of  the  anomalous  heat 
flux  -Pq<s^v^j,>  on  the  parameter  D,  which  we  have  been  after  in  this 
section. 

In  summary,  we  have  derived  the  dependence  of  the  anomalous  heat 
flux  on  the  parameter  D,  which  is  proportional  to  the  mean  pressure 
gradient  p^.  Since  this  derivation  is  based  on  a  small  amplitude 
expansion,  the  result  is  correct  only  when  D  is  slightly  larger  than 
the  linear  stability  limit  D^.  In  this  region,  it  is  found  that  the 
anomalous  heat  flux  is  linear  with  (D-D^^).  The  dependence  on  other 
parameters  R,  M  and  x  af"e,  however,  not  apparent  in  our  results  since 
the  expression  of  the  anomalous  heat  transport  given  in  Eq.  (7-20)  is 
a  functional  of  the  linear  marginal  modes,  which  depend  on  R,  M  and  x 
implicitly. 

When  the  parameter  D  is  much  larger  than  the  linear  stability 
limit  D^^,  the  nonlinearity  of  the  system  is  so  strong  that  we  need 
numerical  calculations  to  estimate  the  anomalous  heat  transport.  In 
the  next  section,  we  will  show  some  results  of  the  numerical 
simulations  and  discuss  the  validity  of  the  analytical  results  of  Eq. 
(7-20). 
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VIII.   NUMERICAL  CALCULATIONS  OF  THE  REDUCED  EQUATIONS 

In  this  section,  the  results  of  numerical  calculations  of  the 
reduced  equations  are  presented  and  some  properties  of  the  anomalous 
heat  transport  are  discussed.  Here  we  only  deal  with  the  system  (T), 
i.e.,  a  tokamak  plasma,  in  which  the  anomalous  heat  transport  is 
decoupled  from  the  dynamo  effect.  It  is  shown  that  steady  convection 
within  a  boundary  layer  is  attained  and  the  anomalous  heat  transport 
arising  from  this  convection  weakly  depends  on  the  collisional  heat 
conductivity.  Therefore,  when  the  effect  of  the  collisional  heat 
conduction  is  small,  the  total  heat  transport  is  dominated  by  the 
anomalous  one  caused  by  the  convective  cells.  It  is  also  shown  that 
the  analytical  estimation  of  the  anomalous  heat  transport  given  by  Eq. 
(7-20)  agrees  well  with  the  results  of  the  numerical  calculations  near 
the  linear  stability  limit  D^. 

The  nonlinear  initial  value  code  to  solve  the  reduced  equations 
(T)  is  developed  from  the  HIB  code  of  Strauss.^'  The  boundary 
conditions  employed  in  our  numerical  code  are  (I),  (II')  and  (III) 
given  in  Section  IV;  the  equations  are  solved  in  the  domain  of  the  slab 
geometry,  i.e.  |x|  <  6^,  |y|  <  6y  and  |z|  <  L^,  which  corresponds  to  a 
boundary  layer  in  the  plasma.  The  numerical  method  used  in  the  code  is 
a  finite  difference  scheme  for  the  variable  x,  Fourier  component 
representations  for  the  variables  y  and  z,  and  the  predictor-corrector 
method  for  time  t.  Since  we  consider  the  fluctuations  on  a  single 
rational  surface,  the  z-dependence  of  the  solutions  is  ignored,  and 
nonlinear  interactions  of  the  modes  with  the  same  hellcity  are 
calculated.  Moreover,  for  simplicity,  only  the  Fourier  cosine 
components  for  (|)  and  p  and  the  Fourier  sine  components  of  A  are  taken 
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into  consideration:   <ti  •  l^  <t>^i^)   cos  (my/yQ),   P  -  Im  Pj^Cx)    cos  (my/yQ), 
and  ^  '  l^  \^^^   2^"   (my/yQ),   where  y^  -  6   /it.     Therefore,     as     in     the 
previous     section,      the  modes  considered  here  are  a  nonlinear  extension 
of   the  normal  mode  representations  of   the  linear  resistive  g-mode. 

Schematically,    the  numerical  algorithm  of  the  system  (T)     employed 
in  our  code   is  summarized  as   follows: 


A*    >   A"  *   hAT[-[^".A"L   *   "x  JH.  ;"  *  I  A, a!) 
m         m  *•    1-^    '      -"m  y     ^n       r         m'' 

p:-;s*hAx(-[;n.-p"], -Df  ;s*xA,p;] 


Yo  '^     r 


^m  m 


*  A  —  A  . 

Here  At  denotes  the  time  step  size,  a"  and  A^  denote  A  and  Ajjj(x)  at 
n-th  time  step,  respectively,  so  that  a"  -  I„  A^  sinCmy/y^),  [  ,  J^^ 
indicates  the  m-th  Fourier  component  of  |  ,  |  and  A  -  d  /dx  - 
(m/y  )^.  The  intermediate  time  step  for  the  predictor-corrector 
advance  is  chosen  by  taking  h  to  be  0.65.  Since  the  implicit  time- 
difference  scheme  is  employed  for  the  operator  A  in  the  right-hand- 
sides  of  these  equations,  operators  of  the  form  A  +  c,  where  c  is  some 
constant,    are  inverted  at  every    time  step     to     solve     these     equations. 
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This     implicit     scheme     allows     us  to  avoid   the  Courant-Friedrichs-Lewy 

22 
condition       and   to   take  a  larger  step  size  At. 

We  now  take  the  small  scale  parameter  6  introduced  in  Section  III 
to  be  1//r,  where  R  is  the  magnetic  Reynolds  number.  In  terms  of  this 
6,  the  following  size  of  the  domain  is  used  for  our  boundary 
conditions:  6^/&  -  50  and  6y/6  -  2ti  (therefore  y^  -  26  in  the  equations 
above).  As  we  will  see  later,  the  solutions  of  these  equations  decay 
rapidly  as  |x|  -►  6^.  It  is  also  confirmed  numerically  that  these 
solutions  hardly  depend  on  the  choice  of  &^  if  6^/6  is  taken  to  be  50 
or  larger.  Therefore  we  use  the  averaging  operator  <  >  defined  in  Eq. 
(^-3)  with  A  -  6  in  order  to  estimate  various  averaged  quantities 
numerically.  In  fact,  it  is  observed  in  these  numerical  computations 
that  the  typical  mode  width  in  the  x-direction  A  is  of  the  same  order 
of  the  one  in   the  y-direction  &    . 

The  boundary  conditions  (I)  and  (II')  imposed  here  depend  on  the 
magnetic  Reynolds  number  R  since  6„  and  6„  are  proportional  to 
6  -  1//r.  However  after  the  scale  transformation  (3-23)  is  applied  to 
this  system,    this  dependence  disappears.      As     we     discussed     after     Eq. 

(3-24)      in     Section     III,      therefore,    the  solutions  of  this  transformed 

~      2      ~      2  ~ 

system   i.e.      A/6    ,    (t>/6     and  p/6,  depend  only  on  the     parameters     D,      M^ 

(-   y  /hPq)   and  Kp   (•    (Y-1 )<   /hPq)   but  not  on  R.      Here  we  call  M^  and  Kp 

the     normalized     viscosity     and     the     normalized       heat       conductivity, 

respectively.       Therefore,      In     actual  numerical  calculations,   we  solve 

this   transformed  system  and  leave   the  magnetic   Reynolds  number  R  as     an 

undetermined  parameter. 

Typically     we     use     the     following   values  of  the  parameters   in  our 

calculations  unless  otherwise  specified:   D   -  0.20  and  Mp,Kf^  -    1    -    10"^. 
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Here  we  note  that  D  <  0.25  is  the  ideal  stability  condition  by 
Suydam.  ^  We  also  note  that  Mpj  -  Kp  -  6  /m~7m^,  where  6  is  the  plasma 
beta,  m^  and  m^  are  the  masses  of  the  ion  and  the  electron, 
respectively,  if  we  take  y  and  k  in  Eqs.  (2-1)  to  be  the  ion 
viscosity  and  the  ion  heat  conductivity,  respectively,  in  a  collisional 
plasma  given  in  Ref.  12.  As  before,  the  parallel  heat  conductivity 
X  is  taken  to  be  smaller  than  0(1)  so  that  it  does  not  enter  the 
equations.  However  we  point  out  that  some  preliminary  calculations 
including  x,,  show  that  x,,  ^^^  s  stabilizing  effect  and  thus  reduces 
somewhat  the  anomalous  heat  transport  across  the  magnetic  surface.  In 
most  runs  presented  here,  150  mesh  ponts  in  the  x-direction  and  7  modes 
(0  _<  m  <_  5)  in  the  y-direction  are  employed.  It  is  confirmed  by 
varying  the  numbers  of  the  mesh  points  and  the  modes  that  this  is  a 
sufficient  numerical  resolution  to  obtain  the  correct  saturation  levels 
of  the  convection.  The  initial  condition  given  to  the  calculations  is 
a  sufficiently  small  perturbation  of  the  m  -  1  component  of  ip. 

Figures  1(a)-1(e)  show  the  Fourier  components  of  A,  A  A,  (p ,  A  0 
and  p,  respectively,  of  the  typical  steady  convection  of  a  tokamak 
obtained  from  Eqs.  (T).  The  dominant  mode  is  the  m  -  1  mode  for  all 
the  quantities  except  for  p,  whose  dominant  mode  is  the  m  -  0  mode.  It 
is  seen,  however,  that  there  are  many  large  modes  other  than  the  m  -  1 
mode,  particularly  in  the  current  density  A  A,  the  vorticity  A  $  and 
the  pressure  p.  Therefore  the  nonlinearity  of  the  system  plays  an 
important  role  in  the  state  of  the  steady  convection.  The  dominance  of 
the  m  -  0  mode  of  p  shows  the  flattening  of  the  mean  pressure  gradient, 
which  is  what  one  expects  a  pressure-driven  instability  to  do.  Figure 
2  shows  the  level  curves  of  $,   i.e.   the  stream  lines  of  the 


-68- 
convection,   corresponding   to  Fig.      1(c).      It  Is  seen  from  these  figures 
that,    as  we  have  expected,    the  modes  decay  rapidly  as    |x|    -►  •     and     the 
typical     widths     of     the  mode   in   the  x-direction  and   in   the  y-dlrection 
are  of   the  same  order. 

The  anomalous  heat  transport  across  the  magnetic  surface  is     often 
described  by   the  Nusselt  number, 


This      is     the     ratio     of     the     actual     heat     flux,    i.e.      the  sum  of  the 

I 

collisional  heat  flux  <  Jq  and  the  anomalous  heat  flux     "Po'^^i^ir''     '^^^ 

to  the  convection,  to  the  collisional  heat  flux  <.Tq.  We  note  that  the 
Nusselt  number  Nu  is  invariant  under  the  scale  transformation  (2-23), 
so  that  Nu  also  depends  only  on  D,  Mp  and  Kp.  Figure  3  shows  the  time 
evolution  of  the  Nusselt  number  of  the  modes  whose  final  state  is 
described  in  Figs.  1(a) -1(e).  It  is  seen  that  the  modes  grow  linearly 
at  the  initial  stage  and  reach  the  steady  convection  state  after  about 
IOOt^,  where  t^  is  the  Alfven  time  defined  as  rQ/p^/(Bg | a | )  (see  Eqs. 
(3-l6)j.  The  curves  (b)  in  Figs.  4  and  5  show  the  dependence  of  the 
Nusselt  number  on  the  parameter  D  obtained  from  the  simulation.  In 
Fig.  i<,  the  normalized  heat  conductivity  Kp  -  10  ,  and  in  Fig.  5, 
Kp  -  10"  .  The  two  straight  lines  (a)  in  these  figures  describe  our 
analytical     results,        calculated        from       Eq.  (7-20).  In       these 

calculations,  the  constant  D2  in  Eq.  (7-20)  was  obtained  through  Eq. 
(7-19)  by  solving  Eq.  (7-11)  numerically.  The  linear  marginal  mode 
solutions  appearing    in  the  right-hand  side  of  Eq.      (7-11)    were  obtained 
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from  the  linear  version  of  the    initial  value  code   for  Eqs.      (T).      It   is 
seen     that  these  Nusselt  numbers  obtained   from  the  bifurcation  analysis 
are   in  good  agreement  with   the  results  of   the  nonlinear  simulation  near 
the  critical  point   D^. 

We  now  seek  the  dependence  of  the  Nusselt  number  Nu  on  the 
normalized  diffusion  coefficients  Mp  and  Kp  for  a  fixed  D.  Figure  6 
shows  that  the  anomalous  heat  flux  "Po^'^l^lr''  vai'ies  weakly  with 
normalized  heat  conductivity  Kp.  Here  the  nondimensional  quantity  Jg 
denotes  the  anomalous  heat  flux  "Po'^^i^ir'*  measured  by  the  quantity 
tiPqIq,  i.e.,  Jg  -  -Pq<s^  v^|,>/(npQTQ) .  This  weak  dependence  of  J^  on  Kp 
is  the  reason  that  the  Nusselt  number  Nu  increases  with  Kp  almost 
linearly.  In  a  realistic  parameter  range  for  fusion  experiments  such 
that  D  -  0.2  and  Mp^,  Kp  -  1.0  -  10~^,  it  is  seen  that  the  Nusselt 
number  Nu  varies  from  1  to  10,  in  other  words,  the  steady  convection 
enhances  the  heat  transport  up  to  one  order  of  magnitude  over  the 
collisional  heat  conduction.  It  is  also  seen  from  Fig.  6  that  when  Kp 
is  less  than  10~^,  the  anomalous  heat  flux  or  the  saturation  level  of 
the  modes  increases  as  some  power  of  Kp  .  This  corresponds  to  the  fact 
that  there  is  no  saturated  mode  nor  stationary  turbulence  when  Kp  -  0, 
which  was  stated  and  proved  as  Proposition  2  in  Section  VI.  Figure  7 
shows  that  the  anomalous  heat  flux  or  the  Nusselt  number  Nu  with  fixed 
Kp  also  depends  weakly  on  the  normalized  viscosity  Mp. 
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IX.   DISCUSSIONS  AND  CONCLUSIONS 

In  this  thesis,  we  have  considered  the  anomalous  heat  transport 
caused  by  the  resistive  g-mode  fluctuations  in  a  plasma.  The  resistive 
g-mode  is  the  instability  caused  by  pressure  or  temperature  gradients 
acting  against  the  curvature  of  the  magnetic  field  lines  in  a  plasma 
with  finite  resistivity,  analogous  to  the  thermal  instability  causing 
the  Benard  convection  in  a  fluid  heated  form  below  in  a  gravity  field. 
We  have  derived  two  systems  of  nonlinear  equations  describing  such 
fluctuations,  one  for  an  RFP  plasma  and  the  other  for  a  tokamak  plasma, 
under  the  assumption  that  the  fluctuating  quantities  have  much  smaller 
scales  in  space  and  time  than  the  corresponding  mean  quantities.  In 
these  equations,  we  take  into  account  the  effects  of  all  the 
collisional  diffusion  coefficients,  i.e.  resistivity,  viscosity  and 
heat  conductivity.  These  systems  are  an  extension  of  the  reduced 
equations  in  Ref .   1 . 

The  inclusion  of  all  the  collisional  diffusion  coefficients  in  the 
reduced  equations  allows  us  to  derive  some  relations  between  different 
kinds  of  anomalous  transport.  It  is  shown  based  on  these  reduced 
equations  that  the  anomalous  heat  transport  <Siy^>  is  related  to  the 
dynamo   effect  e   in  the  RFP  plasma  through  Eq.    (i<-10),   i.e., 

-*-o  *  '^o*'2l-1'**'^o  '  ~  '  where  A  is  the  rate  of  the  energy 
dissipation  of  the  fluctuations.  Therefore  one  should  expect  that 
there  is  a  large  anomalous  heat  transport  when  there  is  a  strong  dynamo 
activity  in  the  RFP.  On  the  other  hand,  it  is  also  shown  that,  in  the 
tokamak,  the  dynamo  effect  is  smaller  than  the  anomalous  heat  transport 

by  an  order  of  the  inverse  aspect  ratio.   In  this  case,  the  anomalous 

?  2 

heat  flux  can  be  written  as  -K  7T  by  using  a  non-negative  function  K  , 
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which  may  depend  on  the  mean  quantities,  including  7Tq.   In  other 
words,  the  anomalous  transport  behaves  like  the  collisional  heat 
conduction  in  the  sense  that  they  both  transport  heat  in  the  direction 
opposite  to  the  temperature  gradient. 

For  the  resistive  g-mode,  the  free-energy  source  is  the  mean 
pressure  (or  temperature)  gradient  and  the  energy  sink  is  the 
collisional  diffusion.  Therefore  when  the  mean  pressure  gradient 
exceeds  a  critical  value,  modes  localized  on  a  rational  surface  begin 
to  grow.  However  these  modes  eventually  saturate  with  finite  amplitude 
if  the  mean  pressure  gradient  is  not  too  large  compared  to  the  critical 
value,  and  we  have  steady  convection  within  a  boundary  layer  which  is 
analogous  to  the  Benard  convection  in  fluid  dynamics.  This  similarity 
between  the  resistive  g-mode  and  the  Benard  convection  is  examined 
based  on  the  system  (T)  of  equations  for  a  tokamak  plasma. 

Applying  the  nonlinear  stability  analysis  to  the  system  (T),  we 
have  shown  that  there  exists  a  positive  critical  value  D^  of  the 
parameter  D,  where  D  is  proportional  to  the  mean  pressure  gradient, 
such  that  all  perturbations  with  arbitrary  amplitude  decay 
monotonically  in  time  if  D  <  D^.  In  fact,  this  D^  is  obtained  as  the 
energy  stability  limit,  or  the  necessary  and  sufficient  condition  for 
monotonic  stability.  It  is  also  shown  that  the  linear  stability  limit 
D^  of  the  system  (T)  is  larger  than  or  equal  to  the  energy  stability 
limit  D  .  Since  all  perturbations  with  arbitrarily  small  amplitude 
grow  at  least  initially  if  D  >  D^i  and  these  perturbations  are  expected 
to  saturate  eventually  with  finite  amplitude,  we  see  the  similarity 
between  the  resistive  g-mode  of  the  system  (T)  and  the  Benard 
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convection  in  fluid  dynamics.   Here  we  note  that  the  parameter  D 
corresponds  to  the  Raylelgh  number  °  of  the  Benard  convection. 

By  regarding  the  steady  convection  as  a  bifurcation  from  the  null 
solutions  of  the  reduced  equations,  we  apply  the  nonlinear  bifurcation 
analysis  to  the  system  (T)  in  order  to  derive  the  dependence  of  the 
anomalous  heat  transport  on  the  mean  pressure  gradient  or  the  parameter 
D.  In  this  method,  the  reduced  equations  are  expanded  in  terms  of  the 
small  amplitude,  and  the  complete  algorithm  to  determine  all  the  higher 
order  terms  is  obtained.  It  is  shown  by  this  method  that,  to  the 
lowest  order,  the  anomalous  heat  transport  is  proportional  to  the 
difference  of  the  parameter  D  and  its  linear  stability  limit  D^,  that 
is,  (D-D^^). 

The  validity  of  the  analysis  mentioned  above  is  confirmed  by  using 
the  numerical  calculations  of  the  nonlinear  reduced  equations.  These 
numerical  calculations  also  determine  the  dependence  of  the  anomalous 
heat  transport  on  various  parameters  in  the  equations.  It  is  found 
from  the  numerical  calculations  that:  (1)  The  anomalous  heat  flux 
obtained  from  the  bifurcation  an.  ./sis  in  Section  III  is  in  good 
agreement  with  the  one  obtained  from  the  numerical  simulations  near  the 
linear  stability  limit  D^;  (2)  The  anomalous  heat  flux  varies  weakly 
with  the  normalized  diffusion  coefficients  Kp  and  M^;  and  (3)  The 
Nusselt  number  Nu,  which  is  the  ratio  of  the  total  heat  flux  to  the 
colllsional  heat  flux,  becomes  therefore  significantly  large  when  the 
collisional  heat  conduction  is  small.  In  a  realistic  parameter  range 
for  fusion  experiments  such  that  D  -  0.2  and  Mp.Kp^  -  1.0  -  10~^,  the 
Nusselt  number  varies  from  1   to   10,   in  other  words,  the  steady 
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convection  on  each  rational  surface  enhances  the  heat  transport  up  to 
one  order  of  the  magnitude  over  the  collisional  heat  conduction. 

In  this  thesis,  up  to  Section  VI,  we  develop  the  theory  which  is 
applicable  to  both  fully  developed  turbulence  and  some  coherent  motion 
such  as  steady  convection.  However,  in  Sections  VII  and  VIII,  where  we 
actually  estimate  the  anomalous  heat  transport,  we  only  treat  coherent 
motion  of  the  plasma  or  the  steady  convection  on  each  rational  surface: 
interactions  of  modes  on  different  rational  surfaces  are  ignored  and 
only  the  nonlinear  interaction  of  the  modes  with  the  same  helicity  is 
considered.  This  assumption  of  the  coherent  structure  of  the  modes 
significantly  simplifies  the  physical  model  and  allows  a  mathematically 
rigorous  treatment  of  the  system.  However,  for  fully  developed 
turbulence,  the  anomalous  heat  transport  can  have  some  different 
features.  Recently  Carreras  et  al.  -^  have  studied  the  fully  developed 
turbulence  caused  by  the  resistive  g-mode,  based  on  a  set  of  reduced 
equations  similar  to  the  system  (T).  They  estimate  the  anomalous  heat 
transport  analytically  by  using  renormalization  techniques  of 
turbulence,  and  numerically  by  solving  the  reduced  equations  in  the 
entire  plasma  so  that  they  are  able  to  take  into  account  interactions 
of  modes  on  different  rational  surfaces.  It  appears  that  their 
estimate  of  the  anomalous  heat  transport  of  fully  developed  turbulence 
is  about  one  order  of  magnitude  larger  than  our  result  for  the  steady 
convection.  Thus  our  results  in  Sections  VII  and  VIII  give  good 
estimates  of  the  anomalous  heat  in  a  relatively  quiet  plasma,  where  we 
may  expect  the  steady  convection  to  persist,  not  to  develop  into 
turbulence. 
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APPENDIX   A 
REDUCED  EQUATIONS  WITH  PARALLEL  VISCOSITY 

Taking  into  account  the  effect  of  the  parallel  viscosity  u  as 
well  as  the  perpendicular  viscosity  y  and  the  heat  conductivity  <,  we 
will  derive  the  set  of  nonlinear  reduced  equations  of  the  resistive 
g-mode  fluctuations.  In  addition  to  the  orderings  introduced  in 
Section  III,  we  assume  that  u  -  0(1).  We  also  assume  that 
V   'v^    '   0('5)   as   in  Section  III  so  that 

§o  -L  ^    -L 

which  is  derived  from  Eq.  (2-lg),  is  0(6)  and  the  (i,j)  component  of  the 
rate  of  strain  tensor  o,  i.e., 


(0,)^^  -  — Li  *  _li  (A-2) 

«'  ■'■J    3x.    9x. 


is     0(1).        These  orderings  are  consistent  with   the  assumptions  made  in 
Section  II. 

Corresponding    to  Eq.      (3~1),    the  momentum  equation  of  v^    is  given, 
up   to  0(6  ) ,   by 

-   -V(p,    ♦   Bq-B,)    +    (Bq-7)B,    +    (Bi-7)Bq   *    (B,-V)bi    -    (V.n)^    .      (A-3) 
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Here 


(V.n)^  -  -3u,,  Bq(B  *B,  ).V[_1J  *  u,  VA,  -  M  [Av,  +  1   V(V.v,)] 
II         1     B^     II       1       3 


Therefore,  the  perpendicular  component  of  Eq.   (A-3)  becomes 


V  (p^  *  ^o'h    ~   3^11^1^  "  °  ^^"^^ 


to  0(1).   Since  V  'B^  -  0  holds  to  0(1),  we  are  able  to  introduce  two 
scalar  functions  A  and  b, .  such  that 


bi-7Axb-b  b  (A-5) 

.1    ^         -    ||. 


as  in  Section  III.  From  Eq.  (A-4),  we  write  the  parallel  component  of 
the  fluctuating  magnetic  field,  i.e.,  b   -  "Bq'B^  as 

b   -  Pi  ♦  3u  X,  ,  (A-6) 

II    '     ir 

assuming  the  mode  vanishes  away  from  the  rational  surface.  With  the 
scalar  functions  <(>  and  v  introduced  in  Eq.  (3-6),  the  operation 
Bq«V    X    (Eq.    (A-3)]  yields 

Po^  Y  -    (§o*Bi)-7(A  A)    -  2bx7(pQ  *  I  !o^'^|| 

-2b.V[I  b2](a  A)    -   3m,,(Jo-Bo)(Bo*§1|)-''^-T^    *   ^^1*    '  ^*"^^ 

20_L  II  1  B^  1  i 
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which     corresponds  to  Eq.      (3"9).     The  parallel  component  of  Eq.      (A-3) 
becomes 


0^  ^  ■  §o-^Pl    *  ?1|-^^,    *  h'^Po  -   3u     b2(b^.b      ).7(^j 


'1        II        -'      ^°  ||-°-°   -'1  b2' 


+    u    .Bq'^'^I    *    l^.^.^  ^A-8) 


11 


The  variable  p^  in  Eq.    (3-8)  should  be  replaced  by  b    in  our 
notation: 


-k^ll'   ^V?1_^)-^^  -  Yl^-^(Po*?o^  -  §o(^-yi)  -  ^Yii  -  0  .   (A-9) 


The  effect  of  viscosity  does  not  enter  the  pressure  equation  to  its 
lowest  order  under  the  assumption  of  our  orderings.  Therefore  we  may 
use  Eq.   (3-^)  as  the  pressure  equation  here. 

Now  we  express  X ^  in  a  more  convenient  form.   From  Eq.   (A-1),  we 
easily  obtain 


^1  '  "  -^  U^B,  ).7v  -K  v^-VIpq  *   I  B^)]  -  ^  V  .y,  .        (A-10) 
B^        1  ^  ^     1 


Therefore  Eqs.  (A-6)-(A-10 ) ,  along  with  Eqs.  (3-3),  (3-^),  (3-7)  and 
(3-11),  form  the  closed  set  of  equations  for  the  resistive  g-mode 
fluctuations. 

We  now  derive  the  nonlinear  reduced  equations  of  fluctuations  in 
terms  of  the  scaled  variables,  as  in  Section  III.  Let  us  assume  that 
the  plasma   is  confined   in  a  cylinder  and     all     mean     quantities     depend 
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only  on  its  radius  r,  as  before.   In  addition  to  the  scaled  variables 
(3-16),  (3-20),  and  the  parameters  (3-17),  (3-19),  we  use  the  following 
scaled  variables 


2r  /^ 


o  '^0  ,         -       2 


Belop^^'^  •     '"pol^P^^'^ 


2  ^  -  ^^o*^ 

b2|o|3/2      II    '  °"   bJo|5/2       1 


and  the  following  parameters 


-r^Po  -2r  (p  *b2)' 


Pol°l'  B2|a|2 


^               3u  Bgloj 
M  -  =L ,      M   -  U __ 

-L   ^0^6  ^         '  I    ^0^'  ^Po 


1    ""^0 
J  -  - (J^-Bo)  . 


Here  we  note  that  the  parameters  N  and  H  are  not  independent  of  the 
parameters  defined  in  Eqs.   (3*19): 

H  -  S  -  D 

D  -Hie  *1n  . 

Y      2 

It  follows  that  Eqs.   (A-6)  and  (3-11)  become 

b  -  p  -  M  A  (A-1 1 ) 
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and  I  ;  -  p  -  111  f  ,  (A-12) 

respectively.   From  Eq.   (A-10),  we  have 


I   .h.   {l,l\   *||i.    1   a   .  (A-13) 

2  3y       3 

36 
Rewriting  Eq  3.      (A-6)-(A-12) ,    (3-3).    (.3-^),    (3-7)   and    (3-11)      in     terms 
of     the     scaled     variable,      we     obtain     the   following   set  of  equations, 
which,    together  with  Eqs.      (A-1 1 ) -(A-1  3) ,   describe   the  resistive  g-mode 

fluctuations  with  the  parallel  viscosity: 


£^.ii*lz  A 

dx  R     _[ 

36 

±11    -  AZ   A    .    lA,A,A}    -^   *   JM,     flA   *    U.A-}]    *M    Z2    5 

^^1-1  1        5y        II    -  11 

36  36 

d    -  ^,   30        - 

dx  3y 


i_b    -   Hil*i^   *    lA."v}    *    ;    -I^b 
dx  3y  ^1 

36 

A  p   -   -oil  -  1!^     *    (Y-1)(K     Z,,    *  K   Z,]t 
dx    ^  3y  2  ^     II     II  11^ 

£^  -!£  *  iZ.b,,|  *  oi^  *  M    (£H  *1  iZ.In  *M  Z  V  . 
dx      _  n         9y       IP    .     3  ^  ■  '^      j_  _L 

36  36 

These  equations  are   the  extension  of  Eqs.      (R)   given  in  Section  III, 
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APPENDIX   B 
THE   BOUNDARY  CONDITIONS  AND   THE    UNIQUENESS   OF   THE  SOLUTIONS 

In  this  Appendix,  we  will  give  a  brief  discussion  of  the  boundary 
conditions  introduced  in  Section  IV  and  we  will  prove  the  uniqueness  of 
the  solutions  of  the  system  (T)  with  these  boundary  conditions, 
assuming  the  existence  of  the  solutions.  The  objective  of  this 
Appendix  is  to  show  that  the  boundary  conditions  are  the  natural  ones 
for  the  nonlinear  reduced  equations  (R)  and  (T)  although  we  do  not 
intend  to  prove  the  well-posedness  of  the  initial-boundary  value 
problems  of   these  systems. 

We  consider  the  solutions  of  the  system  (R)  [or  (T)j  satisfying 
the  following  conditions:  a6C^(I]),  (p  E  C^  iH) ,  p,f,  v  G  C^(fl)  ,  where  Q 
is   the  bounded  doman  defined  by 

n  -    i(x,y,z)  G  R^    :    |x|    <  &^,  |y|    <   6y    .  0   <  z   <  L^}    .      (B-1 ) 

It  is  easy  to  check  that  the  boundary  conditions  (I),  (II)  (or  II')  and 
(III)  are  necessary  and  sufficient  to  derive  "the  energy  equations" 
(4-6)  and  (4-7)  from  the  system  (R).  [For  the  system  (T),  those 
boundary  conditions  for  (p ,  A  4,  A  and  p  are  necessary  and  sufficient  to 
derive  Eqs.  (4-6)  and  (5-3) j.  For  example,  we  need  to  specify  the 
boundary  values  of  t  <t>  as  well  as  41  in  order  to  eliminate  the  surface 
integral   from  the   following   integral   (Green's   identity) 
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-~p-     -  - 
p  (|iA70  dxdy 

^^      1 


where 


Jfi'^       J_  ■'an'^  dn  J_       dn 


n^   -    l(x,y)  G  r2    :    |x|    <   6^    ,      |y|    <   6yl,  (B-2) 


and  d/dn  denotes  differentiation  in  the  direction  of  the  exterior 
normal  to  the  boundary  3Q^.  This  integration  is  used  to  derive  Eq. 
(4-6). 

We  now  prove  the  uniqueness  of  the  solutions  of  the  system  (T). 
Suppose  that  there  exist  solutions  k^.k^E.  C^Cfl),  (^^,<^2^  C  (Ji)  and 
Pl,P2€  C'^(j;),  satisfying  the  boundary  conditions  (I),  (II)  [or  (II')] 
and  (III)  from  t  -  0  to  t  -  T^.,  such  that  A^  -  ^2'  '^■\  '  '^2  ^^^  ^1  '  ^2 
at  T  -  0.  Let  A*  -  A2"A^  ,  <))*  -  '^2~'^\  ^"'^  P*  "  P2"Pi'  then  these 
variables  satisfy   the   following  equation: 


^  e»(T)    -  -<(^«(^i.^  **}>  -   <P*U*.Pi}>   -   (l*D)<p»  ^>  -   F»(t)    ,      (B-3) 

3y 


where 


e»(T)   -  I   (<|V  A»|2>   ^  <|v  4»|2>   +   <|p»|2>), 


and 


F»(t)   -  -  <|a  A»|2>   ♦  M<|A,(|.»|2>   +   x<|V,P*l^> 
^1  1  1 
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Equation  (B-3)    may  be  derived    in  the  way  similar  to  deriving  Eq.      (6-6) 
from  Eqs.      (T).      We  note  that  e*(0)   -   0  from  the  initial  conditions. 
Using  Schwartz's   inequality,   we  obtain 


<p»  Ml>|    <   (<lp»|2><|y«(2>)l/2 
3y 


<I   (<|p»|2>   *   <1V^*I^>) 


<   e»(T)    .  (B-4) 


Similarly 


i<p.ii:!!i>i  <<ip.ii!i>  ||!^||, 

9x     8y  3x  3y 


<    ||!iL|L  e»(T)    .  (B-5) 

3y 


■•2  /  o  ^ 


where   | | f | |,  -.sup    | f(x) |    for  f  €  C(fi) .     Since  p^  G   C^Cn) ,    there  exists 

X  G  « 
a  positive  constant  C,    such  that 


|<P»U».Pi}>|    <  C^    e»(T)    .  (B-6) 


Taking   into  account  the  following  identity 
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9x      3x    3y      3y      3x   3x      3y   3y 


we  obtain 


3x   3x        3y   3y 


Therefore,  there  exist  constants  C2  and  C2  such  that 


<  C2  e»(T)  ,  (B-7) 

the  derivation  of  which  is  similar  to  the  derivation  of  Eqs.   (B-5)  and 
(B-6). 

Thus,  from  Eqs.   (6-10),  (B-3),  (B-H),    (B-6)  and  (B-7),  we  have 

^  e»(T)  <  [C2  *  C^  *    |UD|  -  ^JtMi) 

-  C3  e»(T)  ,  (B-8) 

where  C^  is  a  constant  (which  is  not  necessarily  positive).   It  follows 
that,  since  e»(0)  -  0, 


0    > 
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[    \:       -  C.  e»(t)J   exp   (-C,t)dt 
0  ut  -J  J 

-   f^  -1  te»(t)    exp   (-C,t)Jdt 

•  e*(T)    exp   (-C3T)    , 

for  all  T  _<  Tf.  Therefore,  we  have  e*(T)  -  0,  i.e.,  A»  -())«-  p*  -  0, 
for  all  T  _<  T^,  which  proves  the  uniquensa  of  the  solutions  of  the 
system  (T). 

To  conclude,  we  have  shown  that  the  boundary  conditions  introduced 
in  Section  IV  are  necessary  and  sufficient  to  derive  all  the  results  in 
this  thesis  and  to  guarantee  the  uniqueness  of  the  solutions  of  the 
system  (T).  We  also  note  that  these  boundary  conditions  are  used  to 
solve  the  reduced  equations  numerically  (Section  VIII).  These  facts 
suggest  that  the  boundary  conditions  imposed  on  the  systems  (R)  and  (T) 
are  the  natural  ones  although  we  did  not  prove  the  well-posedness  of 
the  problem. 
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APPENDIX  C 
ON   THE  MAXIMUM  OF  -<pa0/3y>/F 

In  this  Appendix,  we  will  briefly  discuss  the  existence  of  the 
functions  A^,  (j)^  and  p^  that  maximize  the  functional  -<p3(ti/9y>/F(T) 
[Eq,  (6-7)  j.  The  existence  of  these  functions  is  used  to  prove  the 
corollary  of  Theorem  1  in  Section  VI.  For  this  purpose,  we  need  to 
define  the  function  space  H  introduced  in  the  proof  of  Theorem  1  more 
specifically.  The  symbols,  terminology  and  basic  properties  of 
function  spaces  to  be  used   in   this  Appendix   are  found  in   Ref.      24. 

In  this  thesis,  as  discussed  in  Appendix  B,  we  have  considered  the 
solutions  of  the  system  (T)  satisfying  the  following  conditions: 
AeC^(Ii),  (|)  e  C^(«)  and  peC^U),  where  n  is  defined  by  (B-1  ) . 
However,  since  it  is  obvious  that  the  maximum  of  -<p3(t>/9y>/F  is 
attained  when  < | A  A | ^>  -  0,  i.e.,  A^  ^  0,  because  of  the  boundary 
conditions,  we  only  need  to  take  the  supremum  of  the  following 
functional  of  a  pair   of  functions  u  -    (((),p): 


-  <pii> 

y(u)  - !L_ 


m<|a  0p>  *  x<|7,p|  > 


Moreover,  since  y  does  not  involve  any  z-derivatives  of  <>  and  p,  we  may 
think  of  (p  and  p  as  functions  of  only  x  and  y  without  loss  of 
generality.     Therefore,    we  define  the   function  space  H  as 
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H  -  1  (^,p)  :  I  e  C^(n^),  p  e  C^U^),  and  p.l.l.l   -  0  on  3n^}  ,   (C-1  ) 

equipped  with  the  norm 

ll(i.P)|lH  -  Mk"ll?  *  MpIII)^^^  .  (c-2) 

Here  n^  is  defined  by  (B-2)  and 

0<|a|<k 

for  f  ec'^(fl),  where  ||'||  is  the  L^-norm  and  D°'f  is  the  (weak)  partial 
derivative  of  f.  We  here  imposed  the  boundary  conditions  (I)  and  (II) 
in  Section  IV:  in  order  to  use  (II'),  instead  of  (II),  we  need  to 
change  the  boundary  condition  in  (C-1)  accordingly.  With  these 
definitions,    the  original  maximum  problem  is  reduced  to 


sup         y(u)    . 
u  e  H 

I|u|Ih*0 

We  now  define  the  following   two  spaces; 


H   -   |(«,p)    :    ^e  H^(n^),    pe  H^(n^)   and  p,$,A   $   -   0     on  d^^\ 


equipped  with  the  norm   ||*|Ih  ^^6^^"^'^  ^y  ^Q-      (C-2),    and 


H  -   l(^,p)    :   i  e  H3(n2)^   p  e  H^(n^)   and  p,4i,A  i  -   O     on  an^J 


equipped  with  the  norm 

l|(i.P)||- -  (Ih"ll3  *  MpMi)''^ 
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Here  H^(n^)  is  a  Sobolev  space  over  n^.  Clearly  H  Is  the  completion  of 
H  with  respect  to  the  norm  ||'||pj.  We  claim  here  that  ((^^.p  )  6  H  or, 


PROPOSITION  C:  There  exists  Uq<E  H  such  that 


VCUq)  -        sup        y(u)    .  (c-3) 

u  e  H 

I|u|Ih''0 


PROOF:  As  we  have  shown  in  the  proof  of  Theorem  1 ,  there  exists  the 
least  upper  bound  of  |f|.  Therefore  there  exists  a  sequence  |u^| 
(||"nllH  *  0)  in  H  such  that  4'(u^)  -►  the  right-hand-side  of  Eq.  (C-3) 
as  n  ■►  ".  By  forming  a  bounded  sequence  1  li^/I  |  u^|  |  ^j} ,  it  is  easy  to 
check   the  right-hand-side  of  Eq.      (C-3)    is  equal   to 

sup         y(u)    . 
u  e  H 

I|u|Ih-i 

Hence  the  proposition  follows  from  the  fact  that  the  imbedding  H  -►  H  is 
compact,    which   is  an  easy  consequence  of  Rellich's   theorem. 

(Q.E.D.) 

Thus  we  have  defined  the  function  space  H  introduced  in  Section  VI 
more  specifically  and  we  have  shown  that  the  functions  (4>q.Pq)  H  that 
maximize  f  can  be  approximated  by  functions  in  H  as  closely  as  possible 

with  respect  to  the  norm  ||*||-.   Assuming   that  there  exist  smooth 

H 

enough  solutions  of  the  system  (T)  with  these  (less  smooth)  initial 
conditions  (j)^  and  p^  as  well  as  A^  ^  0,  the  discussion  in  this  Appendix 
completes  the  proof  of  the  corollary  of  Theorem  1 . 
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FIGURE  CAPTIONS 

Fig.  1  A  typical  steady  convection  on  a  rational  surface  obtained  from 
Eqs.  (T)  at  time  x  -  300t;^.  The  profiles  of  the  Fourier 
components  of  (a)  A/6^,  (b)  A  A,  (c)  41/6^,  (d)  A  i,  and  (e)  p/6 
are  shown  as  functions  of  x/6,  where  6  -  /r.  Here  12  modes  are 
included,  Mp  -  1 ,  Kp  -  5x10"^,  and  D  -  0.2. 

Fig.  2  Streamlines,  i.e.  contours  of  constant  41,  of  convection  cells, 
corresponding  to  Fig.  1(c). 

Fig.  3  Time  evolution  of  the  Nusselt  number  Nu  under  the  same 
condition  as  in  Fig.  1. 

Fig.  4  The  Nusselt  number  Nu  as  a  function  of  D.   The  straight  line 

(a)  is  the  analytical  prediction  from  Eq.   (7-20)  and  the  curve 

(b)  is  from  the  simulation  of  Eqs.   (T).   Here  Mp  =  1   and 
Kp^  -  10'^. 

Fig.  5  The  Nusselt  number  Nu  as  a  function  of  D.  The  straight  line 

(a)  is  the  analytical  prediction  from  Eq.   (7-20)  and  the  curve 

(b)  is  from  the  simulation  of  Eqs.   (T).  Here  M^   -  1  and 
Kp  -  10"^  . 

Fig.  6  The  dependence  of  the  anomalous  heat  transport  on  the 
normalized  heat  conductivity  K^.  The  curve  (a)  is  the  Nusselt 
number  Nu  and  the  curve  (b)  is  the  anomalous  heat  flux  Jg  = 
-pQ<s^v^^>/nPoTo.  Here  Mr  -  1  and  D  -  0.2. 

Fig.  7  The  Nusselt  number  Nu  as  a  function  of  the  normalized  viscosity 
Mj^  obtained  from  the  simulation  of  Eqs.  (T).  The  curve  (a)  is 
for  Kp  -  10"^  and  the  curve  (b)  is  for  Kp  -  10"\  Here 
D  -  0.2. 
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